Alpha Diversity
2021-02-17
library(knitr)
options(max.print="75")
opts_chunk$set(echo = TRUE, cache = TRUE, warning = FALSE, prompt=FALSE, tidy = TRUE, comment = NA, message = FALSE)
opts_knit$set(width=75)
# Loading library here
library(tidyverse)
#library(ggpubr)
library(MASS)
library(segmented)
library(agricolae)
# set path
source("~/Dropbox/42-JCHWANG-RUTGERS/Projects-Rutgers/Src/utils.R")Initial setup
Setup color index
# Exposure level
Urban_color <- c("#66c2a5", "#fc8d62", "#8da0cb")
names(Urban_color) <- c("Low", "Medium", "High")
print(Urban_color) Low Medium High
"#66c2a5" "#fc8d62" "#8da0cb" # plot(1:3, 1:3, col = Urban_color, pch = 16, cex = 4) Ethnicity
Ethnicity_color <- gg_color_hue(3)
names(Ethnicity_color) <- c("SANEMA", "YEKWANA", "Visitors")
print(Ethnicity_color) SANEMA YEKWANA Visitors
"#F8766D" "#00BA38" "#619CFF" Sanemas <- c("Chajuranha", "Mosenahanha", "Kuyuwininha", "Shianana-Jiyakwanha", "Washudihanha",
"Sudukuma")
Yekwana <- c("Kanarakuni", "Fiyakwanha")Import the data
# Alpha diversities
alpha = readRDS("../data/alpha_imported_all.RDS")# Mapping This mapping files include all villagers and baseline of visitors of
# all body sites
mt_ms = readRDS("../data/Mapping_MS_Urban.Rds")
Mapping_work <- mt_ms %>% mutate(Urban = ifelse(Ethnicity == "SANEMA", "Low", ifelse(Village ==
"Fiyakwanha", "Low", ifelse(SampleGroup == "Visitors", "High", "Medium"))) %>%
factor(levels = c("Low", "Medium", "High")), Age_num = as.numeric(Age), Age_num = ifelse(is.na(Age_num),
as.numeric(gsub(pattern = "_months", replacement = "", Age, ignore.case = T))/12,
Age_num), Age_grp1 = ifelse(Age_num >= 18, "Adults", "Children"))
Mapping_work$Age_num[Mapping_work$Age == "4_Days"] <- 0
Mapping_work$Age_num[Mapping_work$Age == "1_day"] <- 0Set up the Analysis
Work <- merge(alpha, Mapping_work, by = 1) %>% mutate(Age_grp2 = case_when(Age_num <=
3 ~ "Age_0-3", Age_num <= 8 ~ "Age_3-8", Age_num < 18 ~ "Age_8-18", T ~ "Adults"))# Work$Body_Site %>% factor %>% levels
Indices <- colnames(alpha[2:5])
Yrs <- c("2015", "2016")
BodySites <- c("Feces", "Mouth", "Nose", "Right_Arm", "Right_Hand")Impacts on exposure level
dat_all = list()Feces
Initial test (no need to run)
We first check the faith_pd which described the overall complexity of the microbial community.
BB = "Feces"
AA = "faith_pd"
# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers")
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)Start: AIC=606.04
get(AA) ~ 1
Df Sum of Sq RSS AIC
+ Urban 1 1050.37 6177.0 583.23
+ Age_num 1 418.01 6809.4 598.63
<none> 7227.4 606.04
+ Gender 1 88.54 7138.8 606.09
Step: AIC=583.23
get(AA) ~ Urban
Df Sum of Sq RSS AIC
+ Age_num 1 236.50 5940.5 579.06
+ Gender 1 135.24 6041.8 581.73
<none> 6177.0 583.23
Step: AIC=579.06
get(AA) ~ Urban + Age_num
Df Sum of Sq RSS AIC
+ Age_num:Urban 1 177.10 5763.4 576.28
+ Gender 1 173.42 5767.1 576.38
<none> 5940.5 579.06
Step: AIC=576.28
get(AA) ~ Urban + Age_num + Urban:Age_num
Df Sum of Sq RSS AIC
+ Gender 1 133.17 5630.2 574.58
<none> 5763.4 576.28
Step: AIC=574.58
get(AA) ~ Urban + Age_num + Gender + Urban:Age_num
Df Sum of Sq RSS AIC
<none> 5630.2 574.58
+ Urban:Gender 1 66.683 5563.5 574.70
+ Age_num:Gender 1 27.428 5602.8 575.81anova(regforward.out)Analysis of Variance Table
Response: get(AA)
Df Sum Sq Mean Sq F value Pr(>F)
Urban 1 1050.4 1050.37 28.5436 3.263e-07 ***
Age_num 1 236.5 236.50 6.4268 0.01225 *
Gender 1 173.4 173.42 4.7126 0.03148 *
Urban:Age_num 1 136.9 136.86 3.7190 0.05565 .
Residuals 153 5630.2 36.80
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The best model from 2015 includes Urban, Age, Urban:Age interaction and Gender (in decreasing importance).
# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers")
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)Start: AIC=655.24
get(AA) ~ 1
Df Sum of Sq RSS AIC
+ Urban 1 703.20 5886.2 636.70
+ Age_num 1 219.36 6370.0 651.07
<none> 6589.4 655.24
+ Gender 1 10.27 6579.1 656.95
Step: AIC=636.7
get(AA) ~ Urban
Df Sum of Sq RSS AIC
+ Age_num 1 136.141 5750.1 634.44
<none> 5886.2 636.70
+ Gender 1 27.951 5858.3 637.83
Step: AIC=634.44
get(AA) ~ Urban + Age_num
Df Sum of Sq RSS AIC
+ Age_num:Urban 1 170.585 5579.5 630.96
<none> 5750.1 634.44
+ Gender 1 37.584 5712.5 635.24
Step: AIC=630.96
get(AA) ~ Urban + Age_num + Urban:Age_num
Df Sum of Sq RSS AIC
<none> 5579.5 630.96
+ Gender 1 20.061 5559.4 632.30anova(regforward.out)Analysis of Variance Table
Response: get(AA)
Df Sum Sq Mean Sq F value Pr(>F)
Urban 1 703.2 703.20 22.4340 4.428e-06 ***
Age_num 1 136.1 136.14 4.3433 0.03858 *
Urban:Age_num 1 170.6 170.59 5.4421 0.02077 *
Residuals 178 5579.5 31.35
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The best model from 2016 includes Urban, Age, Urban:Age interaction (in decreasing importance), but Gender was not there.
So we will fitting the data with Urban, Age, Urban:Age interaction, and Gender for both years to keep unity.
# Final model
lm_models <- list()
for (YY in Yrs) {
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers")
fit.final <- lm(as.formula(paste0(AA, " ~ Age_num + Urban + Age_num:Urban + Gender")),
data = Work_dat)
print(paste0("--------Below is ANOVA table for year ", YY))
anova(fit.final) %>% print()
print(paste0("--------Below is coefficient for year ", YY))
summary(fit.final) %>% print()
lm_models[[YY]] <- fit.final
}[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: faith_pd
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 418.0 418.01 11.3594 0.0009503 ***
Urban 1 868.9 868.86 23.6110 2.893e-06 ***
Gender 1 173.4 173.42 4.7126 0.0314847 *
Age_num:Urban 1 136.9 136.86 3.7190 0.0556480 .
Residuals 153 5630.2 36.80
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = as.formula(paste0(AA, " ~ Age_num + Urban + Age_num:Urban + Gender")),
data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-14.1948 -4.6478 0.4276 4.8227 14.0085
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 26.79492 1.04178 25.720 < 2e-16 ***
Age_num 0.04259 0.03749 1.136 0.2578
UrbanMedium -7.04666 1.47745 -4.769 4.27e-06 ***
GenderM -1.87206 0.98408 -1.902 0.0590 .
Age_num:UrbanMedium 0.12776 0.06625 1.928 0.0556 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6.066 on 153 degrees of freedom
Multiple R-squared: 0.221, Adjusted R-squared: 0.2006
F-statistic: 10.85 on 4 and 153 DF, p-value: 9.041e-08
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: faith_pd
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 219.4 219.36 6.9839 0.008962 **
Urban 1 620.0 619.99 19.7390 1.561e-05 ***
Gender 1 37.6 37.58 1.1966 0.275488
Age_num:Urban 1 153.1 153.06 4.8732 0.028563 *
Residuals 177 5559.4 31.41
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = as.formula(paste0(AA, " ~ Age_num + Urban + Age_num:Urban + Gender")),
data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-16.7396 -3.4067 -0.0396 3.6630 16.7508
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 24.03558 0.84472 28.454 < 2e-16 ***
Age_num 0.02299 0.02872 0.800 0.4246
UrbanMedium -6.36322 1.38225 -4.604 7.9e-06 ***
GenderM -0.68035 0.85130 -0.799 0.4253
Age_num:UrbanMedium 0.13204 0.05982 2.208 0.0286 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.604 on 177 degrees of freedom
Multiple R-squared: 0.1563, Adjusted R-squared: 0.1372
F-statistic: 8.198 on 4 and 177 DF, p-value: 4.349e-06Next we perform the broken stick regression based on the selected final model
for (YY in Yrs) {
# YY = '2016'
fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = c(3))
summary(fit.seg)
plot(fit.seg$fitted.values, fit.seg$residuals)
plot(fit.seg$model$faith_pd, fit.seg$fitted.values)
plot(fit.seg$model$Age_num, hatvalues(fit.seg))
plot(fit.seg$model$Age_num, cooks.distance(fit.seg))
plot(hatvalues(fit.seg), fit.seg$residuals)
}
It seems from the residual plots that the point with faith_pd smaller than 10 is dragging the curve down and causing non-random residual.
Work %>% filter(Body_Site == BB, faith_pd < 12, SampleGroup == "Villagers") Row.names faith_pd pielou_e shannon observed_otus BarcodeSequence
1 X15.104 11.52192 0.6547718 3.985899 68 ATTTAGGACGAC
2 X15.144 8.21380 0.5725987 3.248024 51 GGTTCCATTAGG
3 X15.149 10.49707 0.6060992 3.455032 52 CGATAGGCCTTA
LinkerPrimerSequence Barcode_location Tube_ID Year SampleGroup Body_Site
1 CCGGACTACHVGGGTWTCTAAT P03_02C 15-104 2015 Villagers Feces
2 CCGGACTACHVGGGTWTCTAAT P03_03D 15-144 2015 Villagers Feces
3 CCGGACTACHVGGGTWTCTAAT P03_03F 15-149 2015 Villagers Feces
Body_Site_Type Village Ethnicity Date House_ID Family_ID
1 Fecal Kanarakuni YEKWANA 10/1/2015-10/3/2015 2015_13 2015_4
2 Fecal Kanarakuni YEKWANA 10/1/2015-10/3/2015 2015_38_39c 2015_21
3 Fecal Kanarakuni YEKWANA 10/1/2015-10/3/2015 2015_6_7 2015_17
Subject_ID Gender Age Urban Age_num Age_grp1 Age_grp2
1 2015_25 M 1 Medium 1.00 Children Age_0-3
2 2015_108 M 15 Medium 15.00 Children Age_8-18
3 2015_143 M 9_months Medium 0.75 Children Age_0-3
[ reached 'max' / getOption("max.print") -- omitted 10 rows ]It does appears that many baby are outliers, so re-do the analysis without baby.
Modeling without babies
dat_p = list()
dat_p[["2015"]] = list()
dat_p[["2016"]] = list()faith_pd
Overview
BB = "Feces"
AA = "faith_pd"
Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >=
1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) +
geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color),
values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8,
panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA,
title = BB)
print(p)
Modeling 2015
# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)Modeling 2016
# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)Final linear model
# Final linear model
lm_models <- list()
for (YY in Yrs) {
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.final <- lm(paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)
print(paste0("--------Below is ANOVA table for year ", YY))
anova(fit.final) %>% print()
print(paste0("--------Below is coefficient for year ", YY))
summary(fit.final) %>% print()
lm_models[[YY]] <- fit.final
}[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: faith_pd
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 251.8 251.81 7.2084 0.008073 **
Urban 1 760.7 760.71 21.7765 6.735e-06 ***
Gender 1 207.1 207.05 5.9272 0.016083 *
Residuals 150 5239.9 34.93
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-15.2172 -4.9553 0.3526 4.6727 12.2899
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 26.70408 0.96592 27.646 < 2e-16 ***
Age_num 0.06583 0.03033 2.170 0.0316 *
UrbanMedium -4.71461 0.98345 -4.794 3.9e-06 ***
GenderM -2.34055 0.96138 -2.435 0.0161 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.91 on 150 degrees of freedom
Multiple R-squared: 0.1888, Adjusted R-squared: 0.1726
F-statistic: 11.64 on 3 and 150 DF, p-value: 6.705e-07
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: faith_pd
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 74.2 74.17 2.6078 0.1081605
Urban 1 436.0 436.01 15.3293 0.0001297 ***
Gender 1 62.2 62.18 2.1861 0.1410845
Residuals 173 4920.7 28.44
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-14.4679 -3.4498 -0.1077 3.4365 16.8511
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 24.21478 0.78318 30.919 < 2e-16 ***
Age_num 0.03226 0.02421 1.333 0.184
UrbanMedium -3.49149 0.87185 -4.005 9.21e-05 ***
GenderM -1.20086 0.81220 -1.479 0.141
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.333 on 173 degrees of freedom
Multiple R-squared: 0.1042, Adjusted R-squared: 0.08866
F-statistic: 6.708 on 3 and 173 DF, p-value: 0.0002618Next we perform the broken stick regression based on the selected final model
stick_models <- list()
for (YY in Yrs) {
# YY = '2016'
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
summary(fit.seg) %>% print()
plot(fit.seg$fitted.values, fit.seg$residuals)
plot(fit.seg$model$faith_pd, fit.seg$fitted.values)
plot(fit.seg$model$Age_num, hatvalues(fit.seg))
plot(fit.seg$model$Age_num, cooks.distance(fit.seg))
plot(hatvalues(fit.seg), fit.seg$residuals)
anova(fit.seg, lm_models[[YY]]) %>% print()
davies.test(lm_models[[YY]], seg.Z = ~Age_num, k = 10) %>% print()
confint.segmented(fit.seg) %>% print()
qqnorm(fit.seg$residuals)
slope(fit.seg) %>% print()
stick_models[[YY]] <- fit.seg
}
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 7.918 1.69
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 20.2457 2.0376 9.936 < 2e-16 ***
Age_num 1.1480 0.4147 2.768 0.00636 **
UrbanMedium -4.3660 0.9399 -4.645 7.45e-06 ***
GenderM -2.2691 0.9107 -2.492 0.01382 *
U1.Age_num -1.1658 0.4157 -2.804 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.578 on 148 degrees of freedom
Multiple R-Squared: 0.287, Adjusted R-squared: 0.2629
Convergence attained in 2 iter. (rel. change 0)
Analysis of Variance Table
Model 1: faith_pd ~ Age_num + Urban + Gender + U1.Age_num + psi1.Age_num
Model 2: faith_pd ~ Age_num + Urban + Gender
Res.Df RSS Df Sum of Sq F Pr(>F)
1 148 4605.7
2 150 5239.9 -2 -634.19 10.19 7.146e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Davies' test for a change in the slope
data: formula = faith_pd ~ Age_num + Urban + Gender , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 7.5556, n.points = 9, p-value = 0.0001387
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 7.91839 4.57902 11.2578
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 1.148000 0.414700 2.76820 0.328470 1.967500
slope2 -0.017832 0.035963 -0.49584 -0.088899 0.053236
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 28 5.359
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 22.55375 0.94600 23.841 < 2e-16 ***
Age_num 0.18250 0.06320 2.888 0.00438 **
UrbanMedium -3.80783 0.85894 -4.433 1.66e-05 ***
GenderM -1.09063 0.79021 -1.380 0.16933
U1.Age_num -0.31478 0.09054 -3.477 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.184 on 171 degrees of freedom
Multiple R-Squared: 0.1634, Adjusted R-squared: 0.1389
Convergence attained in 3 iter. (rel. change 8.2738e-07)
Analysis of Variance Table
Model 1: faith_pd ~ Age_num + Urban + Gender + U1.Age_num + psi1.Age_num
Model 2: faith_pd ~ Age_num + Urban + Gender
Res.Df RSS Df Sum of Sq F Pr(>F)
1 171 4595.5
2 173 4920.7 -2 -325.12 6.0488 0.002896 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Davies' test for a change in the slope
data: formula = faith_pd ~ Age_num + Urban + Gender , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 29.889, n.points = 9, p-value = 0.005568
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 27.9995 17.4216 38.5775
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 0.18250 0.063195 2.8879 0.05776 0.3072500
slope2 -0.13228 0.064494 -2.0510 -0.25959 -0.0049704Creating dataset for plots
dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers",
Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers",
Age_num >= 1)
dat15 <- cbind(dat15, predict(stick_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(stick_models[["2016"]], newdata = dat16, interval = "confidence"))
dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16observed_otu
Overview
BB = "Feces"
AA = "observed_otus"
Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >=
1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) +
geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color),
values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8,
panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA,
title = BB)
print(p)
# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)# Final linear model
lm_models <- list()
for (YY in Yrs) {
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.final <- lm(paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)
print(paste0("--------Below is ANOVA table for year ", YY))
anova(fit.final) %>% print()
print(paste0("--------Below is coefficient for year ", YY))
summary(fit.final) %>% print()
lm_models[[YY]] <- fit.final
}[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: observed_otus
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 41607 41607 6.7746 0.0101742 *
Urban 1 92047 92047 14.9875 0.0001611 ***
Gender 1 33787 33787 5.5013 0.0203122 *
Residuals 150 921238 6142
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-166.815 -60.028 -7.596 56.206 200.962
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 241.7114 12.8076 18.872 < 2e-16 ***
Age_num 0.8844 0.4022 2.199 0.029421 *
UrbanMedium -52.0959 13.0400 -3.995 0.000101 ***
GenderM -29.8986 12.7473 -2.345 0.020312 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 78.37 on 150 degrees of freedom
Multiple R-squared: 0.1538, Adjusted R-squared: 0.1369
F-statistic: 9.088 on 3 and 150 DF, p-value: 1.449e-05
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: observed_otus
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 3124 3124 0.8962 0.345123
Urban 1 32964 32964 9.4556 0.002447 **
Gender 1 3173 3173 0.9100 0.341438
Residuals 173 603117 3486
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-155.109 -46.839 -0.734 43.746 195.489
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 205.3456 8.6706 23.683 < 2e-16 ***
Age_num 0.1914 0.2680 0.714 0.47600
UrbanMedium -30.2261 9.6523 -3.131 0.00204 **
GenderM -8.5779 8.9919 -0.954 0.34144
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 59.04 on 173 degrees of freedom
Multiple R-squared: 0.06112, Adjusted R-squared: 0.04484
F-statistic: 3.754 on 3 and 173 DF, p-value: 0.01204Next we perform the broken stick regression based on the selected final model
stick_models <- list()
for (YY in Yrs) {
# YY = '2015'
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
summary(fit.seg) %>% print()
plot(fit.seg$fitted.values, fit.seg$residuals)
plot(fit.seg$model$observed_otus, fit.seg$fitted.values)
anova(fit.seg, lm_models[[YY]]) %>% print()
davies.test(lm_models[[YY]], seg.Z = ~Age_num, k = 10) %>% print()
confint.segmented(fit.seg) %>% print()
qqnorm(fit.seg$residuals)
slope(fit.seg) %>% print()
stick_models[[YY]] <- fit.seg
}
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 8.333 1.747
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 170.512 24.597 6.932 1.19e-10 ***
Age_num 12.481 4.284 2.913 0.004131 **
UrbanMedium -48.451 12.635 -3.835 0.000186 ***
GenderM -29.088 12.246 -2.375 0.018818 *
U1.Age_num -12.605 4.307 -2.927 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 75.27 on 148 degrees of freedom
Multiple R-Squared: 0.2298, Adjusted R-squared: 0.2038
Convergence attained in 2 iter. (rel. change 0)
Analysis of Variance Table
Model 1: observed_otus ~ Age_num + Urban + Gender + U1.Age_num + psi1.Age_num
Model 2: observed_otus ~ Age_num + Urban + Gender
Res.Df RSS Df Sum of Sq F Pr(>F)
1 148 838484
2 150 921238 -2 -82755 7.3035 0.0009442 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Davies' test for a change in the slope
data: formula = observed_otus ~ Age_num + Urban + Gender , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 7.5556, n.points = 9, p-value = 0.001613
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 8.3333 4.88077 11.7858
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 12.48100 4.28410 2.91330 4.0150 20.94700
slope2 -0.12465 0.52149 -0.23902 -1.1552 0.90588
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 28.004 6.004
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 188.2708 10.4042 18.096 < 2e-16 ***
Age_num 1.7358 0.6676 2.600 0.010131 *
UrbanMedium -33.4771 9.6220 -3.479 0.000638 ***
GenderM -7.4450 8.8148 -0.845 0.399509
U1.Age_num -3.2363 1.0119 -3.198 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 57.67 on 171 degrees of freedom
Multiple R-Squared: 0.1146, Adjusted R-squared: 0.08873
Convergence attained in 5 iter. (rel. change 4.3811e-06)
Analysis of Variance Table
Model 1: observed_otus ~ Age_num + Urban + Gender + U1.Age_num + psi1.Age_num
Model 2: observed_otus ~ Age_num + Urban + Gender
Res.Df RSS Df Sum of Sq F Pr(>F)
1 171 568752
2 173 603117 -2 -34365 5.166 0.006631 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Davies' test for a change in the slope
data: formula = observed_otus ~ Age_num + Urban + Gender , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 29.889, n.points = 9, p-value = 0.01317
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 28.0042 16.1526 39.8559
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 1.7358 0.66756 2.6002 0.41809 3.0535
slope2 -1.5005 0.76088 -1.9721 -3.00250 0.0014Creating data for plots
dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers",
Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers",
Age_num >= 1)
dat15 <- cbind(dat15, predict(stick_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(stick_models[["2016"]], newdata = dat16, interval = "confidence"))
# lm model children
for (YY in Yrs) {
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1, Age_num <= 8)
fit.child <- lm(paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)
print(paste0("--------Below is ANOVA table for year ", YY))
anova(fit.child) %>% print()
print(paste0("--------Below is coefficient for year ", YY))
summary(fit.child) %>% print()
}[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: observed_otus
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 68041 68041 14.3091 0.0003902 ***
Urban 1 10433 10433 2.1940 0.1443576
Gender 1 9950 9950 2.0924 0.1538140
Residuals 54 256776 4755
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-111.68 -47.16 -15.23 32.11 198.99
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 154.237 25.038 6.160 9.43e-08 ***
Age_num 13.557 3.991 3.396 0.00129 **
UrbanMedium -26.009 18.732 -1.388 0.17069
GenderM -26.505 18.324 -1.447 0.15381
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 68.96 on 54 degrees of freedom
Multiple R-squared: 0.2562, Adjusted R-squared: 0.2148
F-statistic: 6.198 on 3 and 54 DF, p-value: 0.00107
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: observed_otus
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 3546 3546.3 1.2597 0.266855
Urban 1 22170 22169.6 7.8753 0.007036 **
Gender 1 489 488.6 0.1736 0.678674
Residuals 52 146385 2815.1
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-108.472 -37.501 -5.733 28.342 110.281
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 181.711 20.126 9.028 3.14e-12 ***
Age_num 3.761 3.575 1.052 0.29760
UrbanMedium -43.276 15.427 -2.805 0.00706 **
GenderM -6.009 14.424 -0.417 0.67867
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 53.06 on 52 degrees of freedom
Multiple R-squared: 0.1518, Adjusted R-squared: 0.1029
F-statistic: 3.103 on 3 and 52 DF, p-value: 0.03442dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16shannon
BB = "Feces"
AA = "shannon"
Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >=
1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) +
geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color),
values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8,
panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA,
title = BB)
print(p)
# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)# Final linear model
lm_models <- list()
for (YY in Yrs) {
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.final <- lm(paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)
print(paste0("--------Below is ANOVA table for year ", YY))
anova(fit.final) %>% print()
print(paste0("--------Below is coefficient for year ", YY))
summary(fit.final) %>% print()
lm_models[[YY]] <- fit.final
}[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: shannon
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 0.778 0.7776 1.1775 0.279607
Urban 1 5.371 5.3707 8.1321 0.004962 **
Gender 1 2.583 2.5833 3.9115 0.049788 *
Residuals 150 99.064 0.6604
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-2.18307 -0.55013 0.05821 0.61292 1.88644
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.822744 0.132813 43.842 < 2e-16 ***
Age_num 0.003427 0.004171 0.822 0.41256
UrbanMedium -0.399803 0.135222 -2.957 0.00361 **
GenderM -0.261435 0.132188 -1.978 0.04979 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8127 on 150 degrees of freedom
Multiple R-squared: 0.081, Adjusted R-squared: 0.06262
F-statistic: 4.407 on 3 and 150 DF, p-value: 0.005302
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: shannon
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 0.304 0.3042 0.5727 0.4502035
Urban 1 7.123 7.1228 13.4092 0.0003325 ***
Gender 1 0.185 0.1845 0.3473 0.5563921
Residuals 173 91.896 0.5312
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-2.23807 -0.46918 -0.03367 0.54295 1.72129
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.552453 0.107028 51.879 < 2e-16 ***
Age_num 0.001486 0.003308 0.449 0.653878
UrbanMedium -0.439992 0.119146 -3.693 0.000297 ***
GenderM -0.065415 0.110994 -0.589 0.556392
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7288 on 173 degrees of freedom
Multiple R-squared: 0.07649, Adjusted R-squared: 0.06048
F-statistic: 4.776 on 3 and 173 DF, p-value: 0.003181Next we perform the broken stick regression based on the selected final model
stick_models <- list()
for (YY in Yrs) {
# YY = '2015'
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
summary(fit.seg) %>% print()
plot(fit.seg$fitted.values, fit.seg$residuals)
plot(fit.seg$model$shannon, fit.seg$fitted.values)
anova(fit.seg, lm_models[[YY]]) %>% print()
davies.test(lm_models[[YY]], seg.Z = ~Age_num, k = 10) %>% print()
confint.segmented(fit.seg) %>% print()
qqnorm(fit.seg$residuals)
slope(fit.seg) %>% print()
stick_models[[YY]] <- fit.seg
}
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 3.571 1.078
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.7605 0.4957 9.605 <2e-16 ***
Age_num 0.3248 0.2050 1.584 0.1152
UrbanMedium -0.3437 0.1352 -2.542 0.0120 *
GenderM -0.2380 0.1309 -1.818 0.0711 .
U1.Age_num -0.3252 0.2050 -1.586 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7994 on 148 degrees of freedom
Multiple R-Squared: 0.1227, Adjusted R-squared: 0.09308
Convergence attained in 2 iter. (rel. change 4.5034e-16)
Analysis of Variance Table
Model 1: shannon ~ Age_num + Urban + Gender + U1.Age_num + psi1.Age_num
Model 2: shannon ~ Age_num + Urban + Gender
Res.Df RSS Df Sum of Sq F Pr(>F)
1 148 94.567
2 150 99.064 -2 -4.4971 3.5191 0.03213 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Davies' test for a change in the slope
data: formula = shannon ~ Age_num + Urban + Gender , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 7.5556, n.points = 9, p-value = 0.1001
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 3.5712 1.44061 5.7018
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 0.32478000 0.2049900 1.584400 -0.080303 0.7298700
slope2 -0.00043922 0.0044672 -0.098321 -0.009267 0.0083885
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 25.999 11.682
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.441128 0.139278 39.067 < 2e-16 ***
Age_num 0.011852 0.009988 1.187 0.237032
UrbanMedium -0.461718 0.120135 -3.843 0.000171 ***
GenderM -0.057888 0.111406 -0.520 0.604003
U1.Age_num -0.019809 0.012968 -1.528 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.728 on 171 degrees of freedom
Multiple R-Squared: 0.08932, Adjusted R-squared: 0.06269
Convergence *not* attained in 1 iter. (rel. change 9.0543e-07)
Analysis of Variance Table
Model 1: shannon ~ Age_num + Urban + Gender + U1.Age_num + psi1.Age_num
Model 2: shannon ~ Age_num + Urban + Gender
Res.Df RSS Df Sum of Sq F Pr(>F)
1 171 90.619
2 173 91.896 -2 -1.2765 1.2044 0.3024
Davies' test for a change in the slope
data: formula = shannon ~ Age_num + Urban + Gender , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 22.667, n.points = 9, p-value = 0.5631
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 25.9988 2.94015 49.0575
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 0.0118520 0.0099883 1.18660 -0.0078642 0.0315680
slope2 -0.0079572 0.0081803 -0.97273 -0.0241050 0.0081901Create data for plots
dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers",
Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers",
Age_num >= 1)
dat15 <- cbind(dat15, predict(stick_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(stick_models[["2016"]], newdata = dat16, interval = "confidence"))
# lm model children
for (YY in Yrs) {
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1, Age_num <= 8)
fit.child <- lm(paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)
print(paste0("--------Below is ANOVA table for year ", YY))
anova(fit.child) %>% print()
print(paste0("--------Below is coefficient for year ", YY))
summary(fit.child) %>% print()
}[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: shannon
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 4.928 4.9281 8.1245 0.006172 **
Urban 1 0.487 0.4867 0.8023 0.374378
Gender 1 0.094 0.0945 0.1558 0.694633
Residuals 54 32.755 0.6066
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-1.79544 -0.55073 0.02717 0.42599 1.58893
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.06029 0.28279 17.894 <2e-16 ***
Age_num 0.11652 0.04508 2.585 0.0125 *
UrbanMedium -0.18395 0.21157 -0.869 0.3884
GenderM -0.08168 0.20695 -0.395 0.6946
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7788 on 54 degrees of freedom
Multiple R-squared: 0.144, Adjusted R-squared: 0.09642
F-statistic: 3.028 on 3 and 54 DF, p-value: 0.03722
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: shannon
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 0.2542 0.25419 0.5781 0.45050
Urban 1 3.0389 3.03886 6.9109 0.01124 *
Gender 1 0.0061 0.00609 0.0138 0.90677
Residuals 52 22.8653 0.43972
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban + Gender"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-1.87772 -0.38368 -0.04977 0.32432 1.43925
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.65450 0.25154 22.479 <2e-16 ***
Age_num -0.03775 0.04468 -0.845 0.4020
UrbanMedium -0.50693 0.19281 -2.629 0.0112 *
GenderM 0.02122 0.18028 0.118 0.9068
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.6631 on 52 degrees of freedom
Multiple R-squared: 0.1261, Adjusted R-squared: 0.07567
F-statistic: 2.501 on 3 and 52 DF, p-value: 0.06955dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16Data for feces
dat_all[[BB]] = dat_pMouth
dat_p = list()
dat_p[["2015"]] = list()
dat_p[["2016"]] = list()Modeling without babies
faith_pd
Overview
BB = "Mouth"
AA = "faith_pd"
Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >=
1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) +
geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color),
values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8,
panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA,
title = BB)
print(p)
# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)For year 2015, only urban and age is significant in determing the faith_pd
# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)In 2016, In addition to Age_num and Urban, the interaction term between age_num and urban is also significant at 0.05. For consistency, only age_num and urban were included.
# Final linear model
lm_models <- list()
for (YY in Yrs) {
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.final <- lm(paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
print(paste0("--------Below is ANOVA table for year ", YY))
anova(fit.final) %>% print()
print(paste0("--------Below is coefficient for year ", YY))
summary(fit.final) %>% print()
lm_models[[YY]] <- fit.final
}[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: faith_pd
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 326.20 326.20 17.846 3.998e-05 ***
Urban 1 219.03 219.03 11.983 0.0006877 ***
Residuals 161 2942.88 18.28
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-6.638 -2.297 -0.327 1.292 35.558
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.01104 0.63355 18.958 < 2e-16 ***
Age_num 0.07478 0.02165 3.454 0.000707 ***
UrbanMedium -2.38385 0.68866 -3.462 0.000688 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.275 on 161 degrees of freedom
Multiple R-squared: 0.1563, Adjusted R-squared: 0.1458
F-statistic: 14.91 on 2 and 161 DF, p-value: 1.142e-06
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: faith_pd
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 578.15 578.15 35.856 1.616e-08 ***
Urban 1 172.54 172.54 10.700 0.00134 **
Residuals 144 2321.87 16.12
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-8.7790 -2.8766 0.1378 2.7790 12.5137
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.07857 0.59478 18.626 < 2e-16 ***
Age_num 0.11555 0.02022 5.714 6.12e-08 ***
UrbanMedium -2.24572 0.68652 -3.271 0.00134 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.015 on 144 degrees of freedom
Multiple R-squared: 0.2443, Adjusted R-squared: 0.2338
F-statistic: 23.28 on 2 and 144 DF, p-value: 1.739e-09Next we perform the broken stick regression based on the selected final model
stick_models <- list()
for (YY in Yrs) {
# YY = '2015'
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
summary(fit.seg) %>% print()
plot(fit.seg$fitted.values, fit.seg$residuals)
plot(fit.seg$model$faith_pd, fit.seg$fitted.values)
qqnorm(fit.seg$residuals)
if (!"segmented" %in% class(fit.seg)) {
print(paste(YY, BB, AA, paste(YY, BB, AA, "same as lm model", sep = "-"),
sep = "-"))
} else {
anova(fit.seg, lm_models[[YY]]) %>% print()
davies.test(lm_models[[YY]], seg.Z = ~Age_num, k = 10) %>% print()
confint.segmented(fit.seg) %>% print()
stick_models[[YY]] <- fit.seg
slope(fit.seg) %>% print()
}
}
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 17.874 9.239
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.89173 0.92518 13.934 < 2e-16 ***
Age_num -0.02790 0.09509 -0.293 0.769591
UrbanMedium -2.38755 0.68826 -3.469 0.000672 ***
U1.Age_num 0.15067 0.10410 1.447 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.268 on 159 degrees of freedom
Multiple R-Squared: 0.1697, Adjusted R-squared: 0.1488
Convergence attained in 5 iter. (rel. change 1.57e-16)
Analysis of Variance Table
Model 1: faith_pd ~ Age_num + Urban + U1.Age_num + psi1.Age_num
Model 2: faith_pd ~ Age_num + Urban
Res.Df RSS Df Sum of Sq F Pr(>F)
1 159 2896.3
2 161 2942.9 -2 -46.535 1.2773 0.2816
Davies' test for a change in the slope
data: formula = faith_pd ~ Age_num + Urban , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 20.667, n.points = 9, p-value = 0.3251
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 17.8742 -0.371911 36.1202
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 -0.02790 0.095089 -0.29341 -0.215700 0.15990
slope2 0.12277 0.043318 2.83430 0.037222 0.20833
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 64.428 6.366
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.95981 0.60095 18.237 < 2e-16 ***
Age_num 0.12518 0.02154 5.811 3.92e-08 ***
UrbanMedium -2.30255 0.68742 -3.350 0.00104 **
U1.Age_num -0.82889 0.94576 -0.876 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.011 on 142 degrees of freedom
Multiple R-Squared: 0.2563, Adjusted R-squared: 0.2354
Convergence *not* attained in 36 iter. (rel. change -0.004225)
Analysis of Variance Table
Model 1: faith_pd ~ Age_num + Urban + U1.Age_num + psi1.Age_num
Model 2: faith_pd ~ Age_num + Urban
Res.Df RSS Df Sum of Sq F Pr(>F)
1 142 2285.1
2 144 2321.9 -2 -36.817 1.144 0.3215
Davies' test for a change in the slope
data: formula = faith_pd ~ Age_num + Urban , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 66, n.points = 9, p-value = 0.3474
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 64.4275 51.844 77.011
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 0.12518 0.021542 5.81090 0.082596 0.16777
slope2 -0.70370 0.945510 -0.74426 -2.572800 1.16540For mouth, broken stick regression did not significantly better than a simple multiple linear regression. So following graph is made using lm.
Data for plots
dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers",
Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers",
Age_num >= 1)
dat15 <- cbind(dat15, predict(lm_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(lm_models[["2016"]], newdata = dat16, interval = "confidence"))
dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16observed_otu
Overview
BB = "Mouth"
AA = "observed_otus"
Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >=
1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) +
geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color),
values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8,
panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA,
title = BB)
print(p)
# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)# Final linear model
lm_models <- list()
for (YY in Yrs) {
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.final <- lm(paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
print(paste0("--------Below is ANOVA table for year ", YY))
anova(fit.final) %>% print()
print(paste0("--------Below is coefficient for year ", YY))
summary(fit.final) %>% print()
lm_models[[YY]] <- fit.final
}[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: observed_otus
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 9449 9449.1 7.0084 0.0089190 **
Urban 1 20971 20970.6 15.5540 0.0001194 ***
Residuals 161 217068 1348.2
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-55.840 -24.481 -4.688 14.206 210.453
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 87.1007 5.4412 16.008 < 2e-16 ***
Age_num 0.3370 0.1859 1.812 0.071811 .
UrbanMedium -23.3257 5.9144 -3.944 0.000119 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 36.72 on 161 degrees of freedom
Multiple R-squared: 0.1229, Adjusted R-squared: 0.112
F-statistic: 11.28 on 2 and 161 DF, p-value: 2.6e-05
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: observed_otus
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 28639 28638.6 18.8704 2.625e-05 ***
Urban 1 10545 10545.4 6.9485 0.009309 **
Residuals 144 218541 1517.6
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-77.785 -28.554 0.636 27.055 114.412
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 78.1266 5.7704 13.539 < 2e-16 ***
Age_num 0.8092 0.1962 4.125 6.25e-05 ***
UrbanMedium -17.5569 6.6604 -2.636 0.00931 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 38.96 on 144 degrees of freedom
Multiple R-squared: 0.152, Adjusted R-squared: 0.1403
F-statistic: 12.91 on 2 and 144 DF, p-value: 6.968e-06Next we perform the broken stick regression based on the selected final model
stick_models <- list()
for (YY in Yrs) {
# YY = '2015'
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
summary(fit.seg) %>% print()
plot(fit.seg$fitted.values, fit.seg$residuals)
plot(fit.seg$model$observed_otus, fit.seg$fitted.values)
qqnorm(fit.seg$residuals)
if (!"segmented" %in% class(fit.seg)) {
print(paste(YY, BB, AA, "same as lm model", sep = "-"))
} else {
anova(fit.seg, lm_models[[YY]]) %>% print()
davies.test(lm_models[[YY]], seg.Z = ~Age_num, k = 10) %>% print()
confint.segmented(fit.seg) %>% print()
stick_models[[YY]] <- fit.seg
slope(fit.seg) %>% print()
}
}
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 23.714 8.224
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 94.6545 6.9218 13.675 < 2e-16 ***
Age_num -0.4631 0.5465 -0.847 0.398024
UrbanMedium -23.0323 5.8751 -3.920 0.000131 ***
U1.Age_num 1.5114 0.7270 2.079 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 36.45 on 159 degrees of freedom
Multiple R-Squared: 0.1463, Adjusted R-squared: 0.1248
Convergence attained in 1 iter. (rel. change 8.3498e-07)
Analysis of Variance Table
Model 1: observed_otus ~ Age_num + Urban + U1.Age_num + psi1.Age_num
Model 2: observed_otus ~ Age_num + Urban
Res.Df RSS Df Sum of Sq F Pr(>F)
1 159 211292
2 161 217068 -2 -5775.9 2.1732 0.1172
Davies' test for a change in the slope
data: formula = observed_otus ~ Age_num + Urban , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 14.111, n.points = 9, p-value = 0.1755
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 23.7139 7.47148 39.9563
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 -0.46312 0.54649 -0.84744 -1.542400 0.6162
slope2 1.04830 0.48321 2.16940 0.093923 2.0026
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 20 11.171
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 85.60221 8.99102 9.521 <2e-16 ***
Age_num -0.02298 0.93448 -0.025 0.9804
UrbanMedium -15.89449 6.74892 -2.355 0.0199 *
U1.Age_num 1.26284 1.01494 1.244 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 38.93 on 142 degrees of freedom
Multiple R-Squared: 0.1651, Adjusted R-squared: 0.1416
Convergence attained in 1 iter. (rel. change 4.0401e-10)
Analysis of Variance Table
Model 1: observed_otus ~ Age_num + Urban + U1.Age_num + psi1.Age_num
Model 2: observed_otus ~ Age_num + Urban
Res.Df RSS Df Sum of Sq F Pr(>F)
1 142 215168
2 144 218541 -2 -3373.4 1.1132 0.3314
Davies' test for a change in the slope
data: formula = observed_otus ~ Age_num + Urban , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 22.667, n.points = 9, p-value = 0.7082
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 20 -2.08335 42.0833
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 -0.022982 0.93448 -0.024593 -1.87030 1.8243
slope2 1.239900 0.38531 3.217800 0.47816 2.0015Again, linear model is sufficiently good
Data for plots
dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers",
Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers",
Age_num >= 1)
dat15 <- cbind(dat15, predict(lm_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(lm_models[["2016"]], newdata = dat16, interval = "confidence"))
dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16shannon
BB = "Mouth"
AA = "shannon"
Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >=
1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) +
geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color),
values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8,
panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA,
title = BB)
print(p)
# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)# Final linear model
lm_models <- list()
for (YY in Yrs) {
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.final <- lm(paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
print(paste0("--------Below is ANOVA table for year ", YY))
anova(fit.final) %>% print()
print(paste0("--------Below is coefficient for year ", YY))
summary(fit.final) %>% print()
lm_models[[YY]] <- fit.final
}[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: shannon
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 1.100 1.1004 1.5312 0.2177
Urban 1 29.613 29.6126 41.2071 1.466e-09 ***
Residuals 161 115.699 0.7186
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-2.71659 -0.53916 0.02071 0.46772 3.11725
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.4484667 0.1256208 35.412 < 2e-16 ***
Age_num -0.0002609 0.0042929 -0.061 0.952
UrbanMedium -0.8765321 0.1365469 -6.419 1.47e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8477 on 161 degrees of freedom
Multiple R-squared: 0.2098, Adjusted R-squared: 0.2
F-statistic: 21.37 on 2 and 161 DF, p-value: 5.877e-09
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: shannon
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 3.651 3.6509 4.0000 0.04738 *
Urban 1 0.109 0.1088 0.1192 0.73044
Residuals 144 131.430 0.9127
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-2.52186 -0.54464 0.02496 0.70104 1.85779
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.782972 0.141509 26.733 <2e-16 ***
Age_num 0.009464 0.004812 1.967 0.0511 .
UrbanMedium -0.056385 0.163336 -0.345 0.7304
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9554 on 144 degrees of freedom
Multiple R-squared: 0.02781, Adjusted R-squared: 0.01431
F-statistic: 2.06 on 2 and 144 DF, p-value: 0.1312Next we perform the broken stick regression based on the selected final model
stick_models <- list()
for (YY in Yrs) {
# YY = '2015'
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
summary(fit.seg) %>% print()
plot(fit.seg$fitted.values, fit.seg$residuals)
plot(fit.seg$model$shannon, fit.seg$fitted.values)
qqnorm(fit.seg$residuals)
if (!"segmented" %in% class(fit.seg)) {
print(paste(YY, BB, AA, "same as lm model", sep = "-"))
} else {
anova(fit.seg, lm_models[[YY]]) %>% print()
davies.test(lm_models[[YY]], seg.Z = ~Age_num, k = 10) %>% print()
confint.segmented(fit.seg) %>% print()
stick_models[[YY]] <- fit.seg
slope(fit.seg) %>% print()
}
}
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-2.71659 -0.53916 0.02071 0.46772 3.11725
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.4484667 0.1256208 35.412 < 2e-16 ***
Age_num -0.0002609 0.0042929 -0.061 0.952
UrbanMedium -0.8765321 0.1365469 -6.419 1.47e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8477 on 161 degrees of freedom
Multiple R-squared: 0.2098, Adjusted R-squared: 0.2
F-statistic: 21.37 on 2 and 161 DF, p-value: 5.877e-09
[1] "2015-Mouth-shannon-same as lm model"
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 3.867 11.041
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.98951 0.91249 4.372 2.36e-05 ***
Age_num -0.04860 0.37329 -0.130 0.897
UrbanMedium -0.04946 0.16576 -0.298 0.766
U1.Age_num 0.05857 0.37336 0.157 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9617 on 142 degrees of freedom
Multiple R-Squared: 0.02855, Adjusted R-squared: 0.00119
Convergence attained in 3 iter. (rel. change -2.2215e-07)
Analysis of Variance Table
Model 1: shannon ~ Age_num + Urban + U1.Age_num + psi1.Age_num
Model 2: shannon ~ Age_num + Urban
Res.Df RSS Df Sum of Sq F Pr(>F)
1 142 131.33
2 144 131.43 -2 -0.10071 0.0544 0.947
Davies' test for a change in the slope
data: formula = shannon ~ Age_num + Urban , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 22.667, n.points = 9, p-value = 0.7103
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 3.86673 -17.9592 25.6927
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 -0.0486030 0.3732900 -0.1302 -0.78653000 0.689320
slope2 0.0099667 0.0051258 1.9444 -0.00016605 0.020099Again broken stick is not sufficiently better than linear
data for plots
dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers",
Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers",
Age_num >= 1)
dat15 <- cbind(dat15, predict(lm_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(lm_models[["2016"]], newdata = dat16, interval = "confidence"))
dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16Data for mouth
dat_all[[BB]] = dat_pNose
dat_p = list()
dat_p[["2015"]] = list()
dat_p[["2016"]] = list()Modeling without babies
faith_pd
BB = "Nose"
AA = "faith_pd"
Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >=
1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) +
geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color),
values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8,
panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA,
title = BB)
print(p)
# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)For year 2015, only urban and age is significant in determing the faith_pd
# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)In 2016, Age_num and Urban as well.
For consistency, only age_num and urban were included.
# Final linear model
lm_models <- list()
for (YY in Yrs) {
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.final <- lm(paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
print(paste0("--------Below is ANOVA table for year ", YY))
anova(fit.final) %>% print()
print(paste0("--------Below is coefficient for year ", YY))
summary(fit.final) %>% print()
lm_models[[YY]] <- fit.final
}[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: faith_pd
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 784.6 784.55 4.5739 0.034276 *
Urban 1 1613.6 1613.59 9.4072 0.002615 **
Residuals 134 22984.6 171.53
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-18.705 -9.132 -3.084 7.123 44.949
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 13.64549 2.04108 6.685 5.67e-10 ***
Age_num 0.18238 0.06913 2.638 0.00932 **
UrbanMedium 7.18087 2.34124 3.067 0.00262 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 13.1 on 134 degrees of freedom
Multiple R-squared: 0.09448, Adjusted R-squared: 0.08096
F-statistic: 6.991 on 2 and 134 DF, p-value: 0.001295
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: faith_pd
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 797.1 797.06 7.1276 0.008523 **
Urban 1 515.7 515.66 4.6112 0.033550 *
Residuals 135 15096.7 111.83
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-17.219 -8.348 -1.427 5.631 30.107
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14.79556 1.66060 8.910 3.08e-15 ***
Age_num 0.16844 0.05563 3.028 0.00295 **
UrbanMedium 4.11374 1.91571 2.147 0.03355 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 10.57 on 135 degrees of freedom
Multiple R-squared: 0.08, Adjusted R-squared: 0.06637
F-statistic: 5.869 on 2 and 135 DF, p-value: 0.003595Next we perform the broken stick regression based on the selected final model
stick_models <- list()
for (YY in Yrs) {
# YY = '2015'
print(paste0("-------Following is for ", YY, "--------------"))
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
summary(fit.seg) %>% print()
plot(fit.seg$fitted.values, fit.seg$residuals)
plot(fit.seg$model$faith_pd, fit.seg$fitted.values)
qqnorm(fit.seg$residuals)
if (!"segmented" %in% class(fit.seg)) {
print(paste(YY, BB, AA, "same as lm model", sep = "-"))
} else {
anova(fit.seg, lm_models[[YY]]) %>% print()
davies.test(lm_models[[YY]], seg.Z = ~Age_num, k = 10) %>% print()
confint.segmented(fit.seg) %>% print()
stick_models[[YY]] <- fit.seg
slope(fit.seg) %>% print()
}
}[1] "-------Following is for 2015--------------"
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 10.503 1.915
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.2593 3.5611 0.354 0.724183
Age_num 2.0584 0.5730 3.592 0.000461 ***
UrbanMedium 7.4982 2.1762 3.446 0.000764 ***
U1.Age_num -2.1586 0.5797 -3.724 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 12.12 on 132 degrees of freedom
Multiple R-Squared: 0.2359, Adjusted R-squared: 0.2127
Convergence attained in 3 iter. (rel. change 0)
Analysis of Variance Table
Model 1: faith_pd ~ Age_num + Urban + U1.Age_num + psi1.Age_num
Model 2: faith_pd ~ Age_num + Urban
Res.Df RSS Df Sum of Sq F Pr(>F)
1 132 19395
2 134 22985 -2 -3589.7 12.216 1.357e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Davies' test for a change in the slope
data: formula = faith_pd ~ Age_num + Urban , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 14.111, n.points = 9, p-value = 7.96e-05
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 10.5028 6.71477 14.2907
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 2.05840 0.573020 3.5922 0.92492 3.191900
slope2 -0.10021 0.096118 -1.0426 -0.29034 0.089921
[1] "-------Following is for 2016--------------"
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 13 2.234
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.3916 2.7305 1.975 0.0504 .
Age_num 1.3496 0.3351 4.028 9.41e-05 ***
UrbanMedium 4.6226 1.7888 2.584 0.0108 *
U1.Age_num -1.4254 0.3456 -4.125 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.854 on 133 degrees of freedom
Multiple R-Squared: 0.213, Adjusted R-squared: 0.1894
Convergence attained in 1 iter. (rel. change -2.1389e-07)
Analysis of Variance Table
Model 1: faith_pd ~ Age_num + Urban + U1.Age_num + psi1.Age_num
Model 2: faith_pd ~ Age_num + Urban
Res.Df RSS Df Sum of Sq F Pr(>F)
1 133 12914
2 135 15097 -2 -2183.2 11.242 3.082e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Davies' test for a change in the slope
data: formula = faith_pd ~ Age_num + Urban , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 15.444, n.points = 9, p-value = 0.0002563
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 12.9999 8.58038 17.4194
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 1.349600 0.335070 4.02800 0.68689 2.012400
slope2 -0.075798 0.087213 -0.86912 -0.24830 0.096705Stick model is better than linear
data for plots
dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers",
Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers",
Age_num >= 1)
dat15 <- cbind(dat15, predict(stick_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(stick_models[["2016"]], newdata = dat16, interval = "confidence"))
# lm model children between 1-8
for (YY in Yrs) {
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1, Age_num <= 8)
fit.child <- lm(paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
print(paste0("--------Below is ANOVA table for year ", YY))
anova(fit.child) %>% print()
print(paste0("--------Below is coefficient for year ", YY))
summary(fit.child) %>% print()
}[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: faith_pd
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 659.2 659.21 3.8258 0.05643 .
Urban 1 300.2 300.19 1.7422 0.19326
Residuals 47 8098.4 172.31
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-14.891 -7.674 -3.182 2.997 41.943
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.1387 4.7677 0.658 0.5135
Age_num 1.8759 0.8572 2.189 0.0336 *
UrbanMedium 5.0475 3.8241 1.320 0.1933
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 13.13 on 47 degrees of freedom
Multiple R-squared: 0.1059, Adjusted R-squared: 0.06787
F-statistic: 2.784 on 2 and 47 DF, p-value: 0.072
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: faith_pd
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 121.3 121.33 1.1391 0.29166
Urban 1 415.6 415.63 3.9020 0.05452 .
Residuals 44 4686.8 106.52
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-13.915 -6.450 -2.914 3.508 31.106
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.4663 4.1339 1.564 0.1249
Age_num 0.9505 0.7221 1.316 0.1949
UrbanMedium 6.0112 3.0431 1.975 0.0545 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 10.32 on 44 degrees of freedom
Multiple R-squared: 0.1028, Adjusted R-squared: 0.06201
F-statistic: 2.521 on 2 and 44 DF, p-value: 0.09197dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16observed_otu
BB = "Nose"
AA = "observed_otus"
Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >=
1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) +
geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color),
values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8,
panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA,
title = BB)
print(p)
# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)# Final linear model
lm_models <- list()
for (YY in Yrs) {
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.final <- lm(paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
print(paste0("--------Below is ANOVA table for year ", YY))
anova(fit.final) %>% print()
print(paste0("--------Below is coefficient for year ", YY))
summary(fit.final) %>% print()
lm_models[[YY]] <- fit.final
}[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: observed_otus
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 19166 19166 1.9738 0.16236
Urban 1 57605 57605 5.9325 0.01618 *
Residuals 134 1301147 9710
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-113.75 -61.06 -27.67 38.30 433.79
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 76.9081 15.3569 5.008 1.7e-06 ***
Age_num 0.9394 0.5202 1.806 0.0732 .
UrbanMedium 42.9052 17.6153 2.436 0.0162 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 98.54 on 134 degrees of freedom
Multiple R-squared: 0.05571, Adjusted R-squared: 0.04162
F-statistic: 3.953 on 2 and 134 DF, p-value: 0.02147
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: observed_otus
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 24477 24477.4 4.6338 0.03313 *
Urban 1 10229 10229.2 1.9365 0.16634
Residuals 135 713122 5282.4
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-90.69 -54.83 -16.74 27.42 270.39
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 76.7036 11.4132 6.721 4.64e-10 ***
Age_num 0.9089 0.3824 2.377 0.0189 *
UrbanMedium 18.3221 13.1665 1.392 0.1663
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 72.68 on 135 degrees of freedom
Multiple R-squared: 0.04641, Adjusted R-squared: 0.03228
F-statistic: 3.285 on 2 and 135 DF, p-value: 0.04045Next we perform the broken stick regression based on the selected final model
stick_models <- list()
for (YY in Yrs) {
# YY = '2015'
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
summary(fit.seg) %>% print()
plot(fit.seg$fitted.values, fit.seg$residuals)
plot(fit.seg$model$observed_otus, fit.seg$fitted.values)
qqnorm(fit.seg$residuals)
if (!"segmented" %in% class(fit.seg)) {
print(paste(YY, BB, AA, "same as lm model", sep = "-"))
} else {
anova(fit.seg, lm_models[[YY]]) %>% print()
davies.test(lm_models[[YY]], seg.Z = ~Age_num, k = 10) %>% print()
confint.segmented(fit.seg) %>% print()
stick_models[[YY]] <- fit.seg
slope(fit.seg) %>% print()
}
}
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 10 1.86
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -14.592 27.044 -0.540 0.590405
Age_num 15.130 4.352 3.477 0.000687 ***
UrbanMedium 45.538 16.526 2.755 0.006690 **
U1.Age_num -16.118 4.402 -3.661 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 92.05 on 132 degrees of freedom
Multiple R-Squared: 0.1882, Adjusted R-squared: 0.1636
Convergence attained in 6 iter. (rel. change 1.8623e-07)
Analysis of Variance Table
Model 1: observed_otus ~ Age_num + Urban + U1.Age_num + psi1.Age_num
Model 2: observed_otus ~ Age_num + Urban
Res.Df RSS Df Sum of Sq F Pr(>F)
1 132 1118559
2 134 1301147 -2 -182588 10.774 4.634e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Davies' test for a change in the slope
data: formula = observed_otus ~ Age_num + Urban , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 7.5556, n.points = 9, p-value = 0.0002517
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 10 6.3204 13.6796
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 15.130 4.35170 3.4769 6.5225 23.7380
slope2 -0.988 0.72994 -1.3535 -2.4319 0.4559
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 13 2.435
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 16.708 18.984 0.880 0.380400
Age_num 8.445 2.330 3.625 0.000411 ***
UrbanMedium 21.568 12.437 1.734 0.085200 .
U1.Age_num -9.094 2.403 -3.785 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 68.51 on 133 degrees of freedom
Multiple R-Squared: 0.1652, Adjusted R-squared: 0.1401
Convergence attained in 4 iter. (rel. change 1.796e-07)
Analysis of Variance Table
Model 1: observed_otus ~ Age_num + Urban + U1.Age_num + psi1.Age_num
Model 2: observed_otus ~ Age_num + Urban
Res.Df RSS Df Sum of Sq F Pr(>F)
1 133 624263
2 135 713122 -2 -88860 9.4658 0.0001434 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Davies' test for a change in the slope
data: formula = observed_otus ~ Age_num + Urban , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 15.444, n.points = 9, p-value = 0.000653
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 13 8.18352 17.8164
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 8.44460 2.32960 3.6249 3.8367 13.05300
slope2 -0.64936 0.60637 -1.0709 -1.8487 0.55002data for plots
dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers",
Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers",
Age_num >= 1)
dat15 <- cbind(dat15, predict(stick_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(stick_models[["2016"]], newdata = dat16, interval = "confidence"))
# lm model children
for (YY in Yrs) {
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1, Age_num <= 8)
fit.child <- lm(paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
print(paste0("--------Below is ANOVA table for year ", YY))
anova(fit.child) %>% print()
print(paste0("--------Below is coefficient for year ", YY))
summary(fit.child) %>% print()
}[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: observed_otus
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 35060 35060 3.9447 0.05287 .
Urban 1 10828 10828 1.2183 0.27531
Residuals 47 417726 8888
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-115.57 -54.57 -29.85 13.78 298.43
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.7517 34.2415 -0.022 0.9826
Age_num 13.3763 6.1561 2.173 0.0349 *
UrbanMedium 30.3149 27.4648 1.104 0.2753
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 94.28 on 47 degrees of freedom
Multiple R-squared: 0.09898, Adjusted R-squared: 0.06064
F-statistic: 2.582 on 2 and 47 DF, p-value: 0.08635
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: observed_otus
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 737 737.0 0.1750 0.6777
Urban 1 5715 5715.4 1.3571 0.2503
Residuals 44 185304 4211.5
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-66.219 -36.428 -23.078 6.143 207.781
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 41.630 25.993 1.602 0.116
Age_num 2.575 4.541 0.567 0.574
UrbanMedium 22.291 19.135 1.165 0.250
Residual standard error: 64.9 on 44 degrees of freedom
Multiple R-squared: 0.03365, Adjusted R-squared: -0.01028
F-statistic: 0.7661 on 2 and 44 DF, p-value: 0.4709dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16shannon
BB = "Nose"
AA = "shannon"
Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >=
1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) +
geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color),
values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8,
panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA,
title = BB)
print(p)
# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)Urban is not significant in both years
# Final linear model
lm_models <- list()
for (YY in Yrs) {
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.final <- lm(paste0(AA, " ~ Age_num"), data = Work_dat)
print(paste0("--------Below is ANOVA table for year ", YY))
anova(fit.final) %>% print()
print(paste0("--------Below is coefficient for year ", YY))
summary(fit.final) %>% print()
lm_models[[YY]] <- fit.final
}[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: shannon
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 11.523 11.5234 7.4168 0.007316 **
Residuals 135 209.750 1.5537
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-2.4027 -0.7720 -0.1378 0.6977 3.9535
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.84986 0.16264 17.522 < 2e-16 ***
Age_num 0.01765 0.00648 2.723 0.00732 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.246 on 135 degrees of freedom
Multiple R-squared: 0.05208, Adjusted R-squared: 0.04506
F-statistic: 7.417 on 1 and 135 DF, p-value: 0.007316
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: shannon
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 7.564 7.5644 5.8903 0.01653 *
Residuals 136 174.653 1.2842
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-2.4861 -0.7870 -0.1429 0.5431 3.9504
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.893257 0.150401 19.237 <2e-16 ***
Age_num 0.014208 0.005854 2.427 0.0165 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.133 on 136 degrees of freedom
Multiple R-squared: 0.04151, Adjusted R-squared: 0.03447
F-statistic: 5.89 on 1 and 136 DF, p-value: 0.01653Next we perform the broken stick regression based on the selected final model
stick_models <- list()
for (YY in Yrs) {
# YY = '2015'
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
summary(fit.seg) %>% print()
plot(fit.seg$fitted.values, fit.seg$residuals)
plot(fit.seg$model$shannon, fit.seg$fitted.values)
qqnorm(fit.seg$residuals)
if (!"segmented" %in% class(fit.seg)) {
print(paste(YY, BB, AA, "same as lm model", sep = "-"))
} else {
anova(fit.seg, lm_models[[YY]]) %>% print()
davies.test(lm_models[[YY]], seg.Z = ~Age_num, k = 10) %>% print()
confint.segmented(fit.seg) %>% print()
stick_models[[YY]] <- fit.seg
slope(fit.seg) %>% print()
}
}
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 12 2.51
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.94561 0.29349 6.629 7.7e-10 ***
Age_num 0.14682 0.04321 3.398 0.000897 ***
U1.Age_num -0.15485 0.04441 -3.487 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.182 on 133 degrees of freedom
Multiple R-Squared: 0.1598, Adjusted R-squared: 0.1409
Convergence attained in 1 iter. (rel. change 9.1184e-09)
Analysis of Variance Table
Model 1: shannon ~ Age_num + U1.Age_num + psi1.Age_num
Model 2: shannon ~ Age_num
Res.Df RSS Df Sum of Sq F Pr(>F)
1 133 185.90
2 135 209.75 -2 -23.847 8.5303 0.0003269 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Davies' test for a change in the slope
data: formula = shannon ~ Age_num , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 14.111, n.points = 9, p-value = 0.0006814
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 12.0005 7.03577 16.9652
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 0.1468200 0.043210 3.39770 0.061349 0.23228
slope2 -0.0080328 0.010259 -0.78297 -0.028325 0.01226
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 14.965 2.962
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.06595 0.27304 7.567 5.5e-12 ***
Age_num 0.11406 0.03512 3.248 0.00147 **
U1.Age_num -0.12816 0.03639 -3.521 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.067 on 134 degrees of freedom
Multiple R-Squared: 0.1627, Adjusted R-squared: 0.1439
Convergence attained in 3 iter. (rel. change 2.7832e-07)
Analysis of Variance Table
Model 1: shannon ~ Age_num + U1.Age_num + psi1.Age_num
Model 2: shannon ~ Age_num
Res.Df RSS Df Sum of Sq F Pr(>F)
1 134 152.58
2 136 174.65 -2 -22.074 9.6931 0.000117 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Davies' test for a change in the slope
data: formula = shannon ~ Age_num , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 15.444, n.points = 9, p-value = 0.0002751
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 14.965 9.1063 20.8238
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 0.114060 0.0351210 3.2476 0.044596 0.1835200
slope2 -0.014103 0.0095434 -1.4777 -0.032978 0.0047725data for plots
dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers",
Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers",
Age_num >= 1)
dat15 <- cbind(dat15, predict(stick_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(stick_models[["2016"]], newdata = dat16, interval = "confidence"))
dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16
# lm model children
for (YY in Yrs) {
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1, Age_num <= 8)
fit.child <- lm(paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
print(paste0("--------Below is ANOVA table for year ", YY))
anova(fit.child) %>% print()
print(paste0("--------Below is coefficient for year ", YY))
summary(fit.child) %>% print()
}[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: shannon
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 3.638 3.6384 2.7025 0.1069
Urban 1 1.025 1.0249 0.7613 0.3874
Residuals 47 63.275 1.3463
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-1.4719 -0.6517 -0.2165 0.4322 4.5073
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.22962 0.42143 5.291 3.13e-06 ***
Age_num 0.10800 0.07577 1.425 0.161
UrbanMedium -0.29494 0.33802 -0.873 0.387
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.16 on 47 degrees of freedom
Multiple R-squared: 0.06864, Adjusted R-squared: 0.02901
F-statistic: 1.732 on 2 and 47 DF, p-value: 0.188
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: shannon
Df Sum Sq Mean Sq F value Pr(>F)
Age_num 1 0.228 0.22845 0.2025 0.6550
Urban 1 0.979 0.97901 0.8676 0.3567
Residuals 44 49.649 1.12840
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Age_num + Urban"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-1.9722 -0.5915 -0.1220 0.2949 2.9743
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.89677 0.42548 6.808 2.18e-08 ***
Age_num -0.04220 0.07432 -0.568 0.573
UrbanMedium -0.29175 0.31321 -0.931 0.357
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.062 on 44 degrees of freedom
Multiple R-squared: 0.02374, Adjusted R-squared: -0.02063
F-statistic: 0.535 on 2 and 44 DF, p-value: 0.5894Data for nose
dat_all[[BB]] = dat_pRight_Arm
dat_p = list()
dat_p[["2015"]] = list()
dat_p[["2016"]] = list()Modeling without babies
faith_pd
BB = "Right_Arm"
AA = "faith_pd"
Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >=
1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) +
geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color),
values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8,
panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA,
title = BB)
print(p)
# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)For year 2015, only urban and gender is significant in determing the faith_pd
# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)In 2016, Only Urban
For consistency, only gender and urban were included.
# Final linear model
lm_models <- list()
for (YY in Yrs) {
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.final <- lm(paste0(AA, " ~ Urban + Gender"), data = Work_dat)
print(paste0("--------Below is ANOVA table for year ", YY))
anova(fit.final) %>% print()
print(paste0("--------Below is coefficient for year ", YY))
summary(fit.final) %>% print()
lm_models[[YY]] <- fit.final
}[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: faith_pd
Df Sum Sq Mean Sq F value Pr(>F)
Urban 1 9453.0 9453.0 53.661 1.083e-11 ***
Gender 1 978.2 978.2 5.553 0.01965 *
Residuals 161 28361.9 176.2
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Urban + Gender"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-31.938 -8.529 -0.744 7.129 43.738
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 41.131 1.681 24.470 < 2e-16 ***
UrbanMedium 15.519 2.093 7.413 6.62e-12 ***
GenderM 4.918 2.087 2.356 0.0197 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 13.27 on 161 degrees of freedom
Multiple R-squared: 0.2689, Adjusted R-squared: 0.2598
F-statistic: 29.61 on 2 and 161 DF, p-value: 1.123e-11
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: faith_pd
Df Sum Sq Mean Sq F value Pr(>F)
Urban 1 1468.4 1468.41 6.3177 0.01396 *
Gender 1 52.2 52.23 0.2247 0.63676
Residuals 80 18594.1 232.43
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Urban + Gender"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-29.123 -11.062 -2.842 11.114 42.245
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 42.860 2.697 15.895 <2e-16 ***
UrbanMedium 8.388 3.347 2.506 0.0142 *
GenderM 1.640 3.460 0.474 0.6368
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 15.25 on 80 degrees of freedom
Multiple R-squared: 0.0756, Adjusted R-squared: 0.05249
F-statistic: 3.271 on 2 and 80 DF, p-value: 0.0431data for plots
dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers",
Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers",
Age_num >= 1)
dat15 <- cbind(dat15, predict(lm_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(lm_models[["2016"]], newdata = dat16, interval = "confidence"))
dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16observed_otu
BB = "Right_Arm"
AA = "observed_otus"
Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >=
1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) +
geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color),
values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8,
panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA,
title = BB)
print(p)
# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)# Final linear model
lm_models <- list()
for (YY in Yrs) {
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.final <- lm(paste0(AA, " ~ Urban + Gender"), data = Work_dat)
print(paste0("--------Below is ANOVA table for year ", YY))
anova(fit.final) %>% print()
print(paste0("--------Below is coefficient for year ", YY))
summary(fit.final) %>% print()
lm_models[[YY]] <- fit.final
}[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: observed_otus
Df Sum Sq Mean Sq F value Pr(>F)
Urban 1 415511 415511 34.9526 1.959e-08 ***
Gender 1 71908 71908 6.0488 0.01497 *
Residuals 161 1913942 11888
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Urban + Gender"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-243.34 -71.27 -10.54 53.40 350.66
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 300.90 13.81 21.791 < 2e-16 ***
UrbanMedium 103.27 17.20 6.005 1.23e-08 ***
GenderM 42.17 17.14 2.459 0.015 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 109 on 161 degrees of freedom
Multiple R-squared: 0.203, Adjusted R-squared: 0.1931
F-statistic: 20.5 on 2 and 161 DF, p-value: 1.171e-08
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: observed_otus
Df Sum Sq Mean Sq F value Pr(>F)
Urban 1 62837 62837 3.2078 0.07707 .
Gender 1 4018 4018 0.2051 0.65185
Residuals 80 1567098 19589
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Urban + Gender"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-235.71 -114.92 -30.32 101.79 454.29
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 326.32 24.75 13.182 <2e-16 ***
UrbanMedium 54.82 30.73 1.784 0.0783 .
GenderM 14.39 31.76 0.453 0.6518
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 140 on 80 degrees of freedom
Multiple R-squared: 0.04092, Adjusted R-squared: 0.01694
F-statistic: 1.706 on 2 and 80 DF, p-value: 0.188data for plots
dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers",
Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers",
Age_num >= 1)
dat15 <- cbind(dat15, predict(lm_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(lm_models[["2016"]], newdata = dat16, interval = "confidence"))
dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16shannon
BB = "Right_Arm"
AA = "shannon"
Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >=
1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) +
geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color),
values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8,
panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA,
title = BB)
print(p)
# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)Urban is not significant in both years
# Final linear model
lm_models <- list()
for (YY in Yrs) {
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.final <- lm(paste0(AA, " ~ Urban + Gender"), data = Work_dat)
print(paste0("--------Below is ANOVA table for year ", YY))
anova(fit.final) %>% print()
print(paste0("--------Below is coefficient for year ", YY))
summary(fit.final) %>% print()
lm_models[[YY]] <- fit.final
}[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: shannon
Df Sum Sq Mean Sq F value Pr(>F)
Urban 1 6.697 6.6969 11.3078 0.0009644 ***
Gender 1 5.370 5.3699 9.0671 0.0030226 **
Residuals 161 95.350 0.5922
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Urban + Gender"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-2.90249 -0.44609 0.04411 0.51664 1.58334
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.36287 0.09746 65.286 < 2e-16 ***
UrbanMedium 0.42235 0.12138 3.480 0.000646 ***
GenderM 0.36439 0.12101 3.011 0.003023 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7696 on 161 degrees of freedom
Multiple R-squared: 0.1123, Adjusted R-squared: 0.1013
F-statistic: 10.19 on 2 and 161 DF, p-value: 6.824e-05
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: shannon
Df Sum Sq Mean Sq F value Pr(>F)
Urban 1 6.884 6.8840 11.6908 0.000991 ***
Gender 1 0.028 0.0284 0.0482 0.826864
Residuals 80 47.108 0.5888
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Urban + Gender"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-2.99333 -0.34828 -0.00307 0.56205 1.52427
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.44814 0.13572 47.509 < 2e-16 ***
UrbanMedium 0.57661 0.16849 3.422 0.000981 ***
GenderM -0.03821 0.17414 -0.219 0.826864
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7674 on 80 degrees of freedom
Multiple R-squared: 0.128, Adjusted R-squared: 0.1062
F-statistic: 5.869 on 2 and 80 DF, p-value: 0.004183data for plots
dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers",
Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers",
Age_num >= 1)
dat15 <- cbind(dat15, predict(lm_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(lm_models[["2016"]], newdata = dat16, interval = "confidence"))
dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16Data for right arm
dat_all[[BB]] = dat_pRight_Hand
dat_p = list()
dat_p[["2015"]] = list()
dat_p[["2016"]] = list()Modeling without babies
faith_pd
BB = "Right_Hand"
AA = "faith_pd"
Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >=
1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) +
geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color),
values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8,
panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA,
title = BB)
print(p)
# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)In 2016, Only Urban
For consistency, only gender and urban were included.
# Final linear model
lm_models <- list()
for (YY in Yrs) {
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.final <- lm(paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)
print(paste0("--------Below is ANOVA table for year ", YY))
anova(fit.final) %>% print()
print(paste0("--------Below is coefficient for year ", YY))
summary(fit.final) %>% print()
lm_models[[YY]] <- fit.final
}[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: faith_pd
Df Sum Sq Mean Sq F value Pr(>F)
Urban 1 17233.1 17233.1 121.850 < 2.2e-16 ***
Gender 1 1705.5 1705.5 12.059 0.0006631 ***
Age_num 1 588.1 588.1 4.158 0.0430834 *
Residuals 160 22628.6 141.4
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-35.668 -7.504 0.730 6.352 33.437
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 39.10624 1.91115 20.462 < 2e-16 ***
UrbanMedium 20.24935 1.89466 10.688 < 2e-16 ***
GenderM 6.61860 1.87575 3.529 0.000546 ***
Age_num -0.11885 0.05829 -2.039 0.043083 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 11.89 on 160 degrees of freedom
Multiple R-squared: 0.4632, Adjusted R-squared: 0.4531
F-statistic: 46.02 on 3 and 160 DF, p-value: < 2.2e-16
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: faith_pd
Df Sum Sq Mean Sq F value Pr(>F)
Urban 1 7975.2 7975.2 36.5302 1.978e-08 ***
Gender 1 799.9 799.9 3.6638 0.05813 .
Age_num 1 545.3 545.3 2.4976 0.11681
Residuals 113 24669.9 218.3
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-34.010 -9.883 -0.960 11.680 29.544
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 37.56772 2.70294 13.899 < 2e-16 ***
UrbanMedium 16.51829 2.74989 6.007 2.35e-08 ***
GenderM 5.43211 2.82161 1.925 0.0567 .
Age_num -0.13625 0.08621 -1.580 0.1168
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 14.78 on 113 degrees of freedom
Multiple R-squared: 0.2742, Adjusted R-squared: 0.2549
F-statistic: 14.23 on 3 and 113 DF, p-value: 6.252e-08Next we perform the broken stick regression based on the selected final model
stick_models <- list()
for (YY in Yrs) {
# YY = '2015'
print(paste0("-------Following is for ", YY, "--------------"))
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
summary(fit.seg) %>% print()
plot(fit.seg$fitted.values, fit.seg$residuals)
plot(fit.seg$model$faith_pd, fit.seg$fitted.values)
qqnorm(fit.seg$residuals)
if (!"segmented" %in% class(fit.seg)) {
print(paste(YY, BB, AA, "same as lm model", sep = "-"))
} else {
anova(fit.seg, lm_models[[YY]]) %>% print()
davies.test(lm_models[[YY]], seg.Z = ~Age_num, k = 10) %>% print()
confint.segmented(fit.seg) %>% print()
stick_models[[YY]] <- fit.seg
slope(fit.seg) %>% print()
}
}[1] "-------Following is for 2015--------------"
Call:
lm(formula = paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-35.668 -7.504 0.730 6.352 33.437
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 39.10624 1.91115 20.462 < 2e-16 ***
UrbanMedium 20.24935 1.89466 10.688 < 2e-16 ***
GenderM 6.61860 1.87575 3.529 0.000546 ***
Age_num -0.11885 0.05829 -2.039 0.043083 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 11.89 on 160 degrees of freedom
Multiple R-squared: 0.4632, Adjusted R-squared: 0.4531
F-statistic: 46.02 on 3 and 160 DF, p-value: < 2.2e-16
[1] "2015-Right_Hand-faith_pd-same as lm model"
[1] "-------Following is for 2016--------------"
Call:
lm(formula = paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-34.010 -9.883 -0.960 11.680 29.544
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 37.56772 2.70294 13.899 < 2e-16 ***
UrbanMedium 16.51829 2.74989 6.007 2.35e-08 ***
GenderM 5.43211 2.82161 1.925 0.0567 .
Age_num -0.13625 0.08621 -1.580 0.1168
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 14.78 on 113 degrees of freedom
Multiple R-squared: 0.2742, Adjusted R-squared: 0.2549
F-statistic: 14.23 on 3 and 113 DF, p-value: 6.252e-08
[1] "2016-Right_Hand-faith_pd-same as lm model"data for plots
dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers",
Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers",
Age_num >= 1)
dat15 <- cbind(dat15, predict(lm_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(lm_models[["2016"]], newdata = dat16, interval = "confidence"))
dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16observed_otu
BB = "Right_Hand"
AA = "observed_otus"
Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >=
1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) +
geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color),
values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8,
panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA,
title = BB)
print(p)
# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)# Final linear model
lm_models <- list()
for (YY in Yrs) {
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.final <- lm(paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)
print(paste0("--------Below is ANOVA table for year ", YY))
anova(fit.final) %>% print()
print(paste0("--------Below is coefficient for year ", YY))
summary(fit.final) %>% print()
lm_models[[YY]] <- fit.final
}[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: observed_otus
Df Sum Sq Mean Sq F value Pr(>F)
Urban 1 938567 938567 92.5543 < 2e-16 ***
Gender 1 76969 76969 7.5901 0.00655 **
Age_num 1 34657 34657 3.4176 0.06635 .
Residuals 160 1622514 10141
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-253.91 -66.40 -0.44 55.81 354.88
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 286.3689 16.1831 17.696 < 2e-16 ***
UrbanMedium 149.1367 16.0434 9.296 < 2e-16 ***
GenderM 44.5666 15.8833 2.806 0.00564 **
Age_num -0.9124 0.4935 -1.849 0.06635 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 100.7 on 160 degrees of freedom
Multiple R-squared: 0.3929, Adjusted R-squared: 0.3815
F-statistic: 34.52 on 3 and 160 DF, p-value: < 2.2e-16
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: observed_otus
Df Sum Sq Mean Sq F value Pr(>F)
Urban 1 534368 534368 28.009 5.994e-07 ***
Gender 1 69292 69292 3.632 0.05922 .
Age_num 1 51988 51988 2.725 0.10157
Residuals 113 2155844 19078
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-282.28 -99.17 -23.80 102.44 335.26
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 281.1263 25.2675 11.126 < 2e-16 ***
UrbanMedium 135.0262 25.7063 5.253 7.13e-07 ***
GenderM 50.5734 26.3768 1.917 0.0577 .
Age_num -1.3304 0.8059 -1.651 0.1016
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 138.1 on 113 degrees of freedom
Multiple R-squared: 0.2332, Adjusted R-squared: 0.2128
F-statistic: 11.46 on 3 and 113 DF, p-value: 1.293e-06Next we perform the broken stick regression based on the selected final model
stick_models <- list()
for (YY in Yrs) {
# YY = '2015'
print(paste0("-------Following is for ", YY, "--------------"))
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
summary(fit.seg) %>% print()
plot(fit.seg$fitted.values, fit.seg$residuals)
plot(fit.seg$model$observed_otus, fit.seg$fitted.values)
qqnorm(fit.seg$residuals)
if (!"segmented" %in% class(fit.seg)) {
print(paste(YY, BB, AA, "same as lm model", sep = "-"))
} else {
anova(fit.seg, lm_models[[YY]]) %>% print()
davies.test(lm_models[[YY]], seg.Z = ~Age_num, k = 10) %>% print()
confint.segmented(fit.seg) %>% print()
stick_models[[YY]] <- fit.seg
slope(fit.seg) %>% print()
}
}[1] "-------Following is for 2015--------------"
Call:
lm(formula = paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-253.91 -66.40 -0.44 55.81 354.88
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 286.3689 16.1831 17.696 < 2e-16 ***
UrbanMedium 149.1367 16.0434 9.296 < 2e-16 ***
GenderM 44.5666 15.8833 2.806 0.00564 **
Age_num -0.9124 0.4935 -1.849 0.06635 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 100.7 on 160 degrees of freedom
Multiple R-squared: 0.3929, Adjusted R-squared: 0.3815
F-statistic: 34.52 on 3 and 160 DF, p-value: < 2.2e-16
[1] "2015-Right_Hand-observed_otus-same as lm model"
[1] "-------Following is for 2016--------------"
Call:
lm(formula = paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-282.28 -99.17 -23.80 102.44 335.26
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 281.1263 25.2675 11.126 < 2e-16 ***
UrbanMedium 135.0262 25.7063 5.253 7.13e-07 ***
GenderM 50.5734 26.3768 1.917 0.0577 .
Age_num -1.3304 0.8059 -1.651 0.1016
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 138.1 on 113 degrees of freedom
Multiple R-squared: 0.2332, Adjusted R-squared: 0.2128
F-statistic: 11.46 on 3 and 113 DF, p-value: 1.293e-06
[1] "2016-Right_Hand-observed_otus-same as lm model"data for plots
dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers",
Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers",
Age_num >= 1)
dat15 <- cbind(dat15, predict(lm_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(lm_models[["2016"]], newdata = dat16, interval = "confidence"))
dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16shannon
BB = "Right_Hand"
AA = "shannon"
Work_dat <- Work %>% filter(Body_Site == BB, SampleGroup == "Villagers", Age_num >=
1)
p <- ggplot(Work_dat, aes(x = Age_num, y = get(AA), color = Urban, shape = Gender)) +
geom_point(size = 1.5) + scale_shape_manual(values = c(16, 21)) + scale_color_manual(breaks = names(Urban_color),
values = Urban_color) + facet_grid(. ~ Year) + theme_jw() + theme(aspect.ratio = 0.8,
panel.background = element_rect(color = "black")) + labs(x = "Age (Years)", y = AA,
title = BB)
print(p)
# Year 2015
YY = "2015"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)# Year 2016
YY = "2016"
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
# Obtain the full linear model
fit.full <- lm(get(AA) ~ Age_num + Urban + Gender + Age_num:Urban + Age_num:Gender +
Urban:Gender, data = Work_dat)
# Obtain the null model with no predictors, this is necessary for forward
# selection use stepAIC
fit.null <- lm(get(AA) ~ 1, data = Work_dat)
# Forward selection, scope parameters provides the full model
regforward.out <- stepAIC(fit.null, direction = "forward", trace = 4, scope = ~Age_num +
Urban + Gender + Age_num:Urban + Age_num:Gender + Urban:Gender)
anova(regforward.out)# Final linear model
lm_models <- list()
for (YY in Yrs) {
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.final <- lm(paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)
print(paste0("--------Below is ANOVA table for year ", YY))
anova(fit.final) %>% print()
print(paste0("--------Below is coefficient for year ", YY))
summary(fit.final) %>% print()
lm_models[[YY]] <- fit.final
}[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: shannon
Df Sum Sq Mean Sq F value Pr(>F)
Urban 1 19.569 19.5685 40.4076 2.058e-09 ***
Gender 1 2.413 2.4132 4.9832 0.02698 *
Age_num 1 2.366 2.3663 4.8863 0.02849 *
Residuals 160 77.484 0.4843
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-2.20037 -0.41410 0.09709 0.48491 1.29455
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.287567 0.111834 56.222 < 2e-16 ***
UrbanMedium 0.665392 0.110869 6.002 1.26e-08 ***
GenderM 0.251744 0.109762 2.294 0.0231 *
Age_num -0.007539 0.003411 -2.210 0.0285 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.6959 on 160 degrees of freedom
Multiple R-squared: 0.2391, Adjusted R-squared: 0.2248
F-statistic: 16.76 on 3 and 160 DF, p-value: 1.621e-09
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: shannon
Df Sum Sq Mean Sq F value Pr(>F)
Urban 1 13.575 13.5753 11.8432 0.0008121 ***
Gender 1 2.522 2.5217 2.1999 0.1408027
Age_num 1 7.293 7.2927 6.3622 0.0130500 *
Residuals 113 129.527 1.1463
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-4.2594 -0.3309 0.1945 0.5960 1.6278
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.205022 0.195854 31.682 < 2e-16 ***
UrbanMedium 0.666162 0.199256 3.343 0.00112 **
GenderM 0.306865 0.204453 1.501 0.13617
Age_num -0.015757 0.006247 -2.522 0.01305 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.071 on 113 degrees of freedom
Multiple R-squared: 0.153, Adjusted R-squared: 0.1305
F-statistic: 6.802 on 3 and 113 DF, p-value: 0.000295Next we perform the broken stick regression based on the selected final model
stick_models <- list()
for (YY in Yrs) {
# YY = '2015'
print(paste0("-------Following is for ", YY, "--------------"))
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup == "Villagers",
Age_num >= 1)
fit.seg <- segmented(lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
summary(fit.seg) %>% print()
plot(fit.seg$fitted.values, fit.seg$residuals)
plot(fit.seg$model$shannon, fit.seg$fitted.values)
qqnorm(fit.seg$residuals)
if (!"segmented" %in% class(fit.seg)) {
print(paste(YY, BB, AA, "same as lm model", sep = "-"))
} else {
anova(fit.seg, lm_models[[YY]]) %>% print()
davies.test(lm_models[[YY]], seg.Z = ~Age_num, k = 10) %>% print()
confint.segmented(fit.seg) %>% print()
stick_models[[YY]] <- fit.seg
slope(fit.seg) %>% print()
}
}[1] "-------Following is for 2015--------------"
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 4.314 1.83
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.7584 0.3410 16.889 < 2e-16 ***
UrbanMedium 0.6897 0.1109 6.217 4.32e-09 ***
GenderM 0.2639 0.1092 2.416 0.0168 *
Age_num 0.1338 0.1195 1.120 0.2643
U1.Age_num -0.1445 0.1195 -1.210 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.6912 on 158 degrees of freedom
Multiple R-Squared: 0.2587, Adjusted R-squared: 0.2352
Convergence attained in 2 iter. (rel. change 1.8799e-16)
Analysis of Variance Table
Model 1: shannon ~ Urban + Gender + Age_num + U1.Age_num + psi1.Age_num
Model 2: shannon ~ Urban + Gender + Age_num
Res.Df RSS Df Sum of Sq F Pr(>F)
1 158 75.492
2 160 77.484 -2 -1.9924 2.085 0.1277
Davies' test for a change in the slope
data: formula = shannon ~ Urban + Gender + Age_num , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 7.5556, n.points = 9, p-value = 0.4543
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 4.31412 0.699437 7.9288
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 0.133820 0.1194600 1.1202 -0.102120 0.3697600
slope2 -0.010718 0.0038693 -2.7701 -0.018361 -0.0030762
[1] "-------Following is for 2016--------------"
***Regression Model with Segmented Relationship(s)***
Call:
segmented.lm(obj = lm_models[[YY]], seg.Z = ~Age_num, psi = 9)
Estimated Break-Point(s):
Est. St.Err
psi1.Age_num 45 8.018
Meaningful coefficients of the linear terms:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.34024 0.21078 30.080 < 2e-16 ***
UrbanMedium 0.69974 0.19811 3.532 0.000602 ***
GenderM 0.27327 0.20319 1.345 0.181383
Age_num -0.02709 0.00933 -2.904 0.004452 **
U1.Age_num 0.07724 0.04678 1.651 NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.06 on 111 degrees of freedom
Multiple R-Squared: 0.184, Adjusted R-squared: 0.1473
Convergence attained in 1 iter. (rel. change 2.109e-07)
Analysis of Variance Table
Model 1: shannon ~ Urban + Gender + Age_num + U1.Age_num + psi1.Age_num
Model 2: shannon ~ Urban + Gender + Age_num
Res.Df RSS Df Sum of Sq F Pr(>F)
1 111 124.78
2 113 129.53 -2 -4.7481 2.1119 0.1258
Davies' test for a change in the slope
data: formula = shannon ~ Urban + Gender + Age_num , method = lm
model = gaussian , link = identity
segmented variable = Age_num
'best' at = 46.889, n.points = 9, p-value = 0.2338
alternative hypothesis: two.sided
Est. CI(95%).low CI(95%).up
psi1.Age_num 44.9996 29.1121 60.8871
$Age_num
Est. St.Err. t value CI(95%).l CI(95%).u
slope1 -0.027091 0.0093302 -2.9036 -0.045580 -0.0086029
slope2 0.050152 0.0458170 1.0946 -0.040638 0.1409400data for plots
dat15 <- Work %>% filter(Body_Site == BB, Year == "2015", SampleGroup == "Villagers",
Age_num >= 1)
dat16 <- Work %>% filter(Body_Site == BB, Year == "2016", SampleGroup == "Villagers",
Age_num >= 1)
dat15 <- cbind(dat15, predict(lm_models[["2015"]], newdata = dat15, interval = "confidence"))
dat16 <- cbind(dat16, predict(lm_models[["2016"]], newdata = dat16, interval = "confidence"))
dat_p[["2015"]][[AA]] = dat15
dat_p[["2016"]][[AA]] = dat16Data for right hand
dat_all[[BB]] = dat_pStats check
Stats on the urban, age, and gender
Indices <- colnames(alpha[c(2, 4:5)])
Yrs <- c("2015", "2016")
BodySites <- c("Feces", "Mouth", "Nose", "Right_Arm", "Right_Hand")
lm_results <- list()
duncan_results <- list()
for (AA in Indices) {
for (BB in BodySites) {
for (YY in Yrs) {
Work_dat <- Work %>% filter(Body_Site == BB, Year == YY, SampleGroup ==
"Villagers", Age_num >= 1)
fit.final <- lm(paste0(AA, " ~ Urban + Gender + Age_num"), data = Work_dat)
print(paste0("--------Below is ANOVA table for year ", YY))
anova(fit.final) %>% print()
print(paste0("--------Below is coefficient for year ", YY))
summary(fit.final) %>% print()
lm_results[[paste(AA, BB, YY, sep = "-")]] <- summary(fit.final)$coefficients[-1,
] %>% as.data.frame() %>% rownames_to_column(var = "Factor") %>%
mutate(Year = YY, Body_Sites = BB, Index = AA)
fit.final.village <- lm(paste0(AA, " ~ Village + Gender + Age_num"),
data = Work_dat)
duncan_results[[paste(AA, BB, YY, sep = "-")]] <- duncan.test(fit.final.village,
trt = "Village", console = T)$groups %>% rownames_to_column(var = "Village") %>%
mutate(Year = YY, Body_Sites = BB, Index = AA)
colnames(duncan_results[[paste(AA, BB, YY, sep = "-")]])[2] <- "Value"
}
}
}
lm_results_final <- do.call("rbind", lm_results)
lm_results_final$fdr <- p.adjust(lm_results_final$`Pr(>|t|)`, method = "fdr")
duncan_results_final <- do.call("rbind", duncan_results)
write.table(lm_results_final, file = "../data/alpha_lm_stats.txt", sep = "\t", row.names = F,
quote = F)
# clip <- pipe('pbcopy', 'w') write.table(lm_results_final, file = clip,
# row.names = F, quote = F, sep = '\t') close(clip)within Sanema village effect
Indices <- colnames(alpha[c(2, 4:5)])
Yrs <- c("2015", "2016")
BodySites <- c("Feces", "Mouth", "Nose", "Right_Arm", "Right_Hand")
lm_results <- list()
anova_results <- list()
for (AA in Indices) {
for (BB in BodySites) {
for (YY in Yrs) {
Work_dat <- Work %>% filter(Ethnicity == "SANEMA", Body_Site == BB, Year ==
YY, SampleGroup == "Villagers", Age_num >= 1)
fit.final <- lm(paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
print(paste0("--------Below is ANOVA table for year ", YY))
anova(fit.final) %>% print()
print(paste0("--------Below is coefficient for year ", YY))
summary(fit.final) %>% print()
lm_results[[paste(AA, BB, YY, sep = "-")]] <- summary(fit.final)$coefficients[-1,
] %>% as.data.frame() %>% rownames_to_column(var = "Factor") %>%
mutate(Year = YY, Body_Sites = BB, Index = AA)
anova_results[[paste(AA, BB, YY, sep = "-")]] <- as.data.frame(anova(fit.final))[1:3,
4:5] %>% rownames_to_column(var = "Factor") %>% mutate(Year = YY,
Body_Sites = BB, Index = AA)
}
}
}[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: faith_pd
Df Sum Sq Mean Sq F value Pr(>F)
Village 4 267.74 66.936 1.9358 0.1119
Gender 1 146.23 146.232 4.2292 0.0428 *
Age_num 1 61.98 61.976 1.7924 0.1842
Residuals 85 2939.03 34.577
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-14.679 -4.596 1.244 4.749 9.781
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 30.05499 2.08522 14.413 <2e-16 ***
VillageKuyuwininha -3.45478 2.24784 -1.537 0.1280
VillageMosenahanha -0.04427 2.55933 -0.017 0.9862
VillageShianana-Jiyakwanha -3.08565 2.18967 -1.409 0.1624
VillageWashudihanha -4.49443 2.27934 -1.972 0.0519 .
GenderM -2.90774 1.28950 -2.255 0.0267 *
Age_num 0.05008 0.03741 1.339 0.1842
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.88 on 85 degrees of freedom
Multiple R-squared: 0.1394, Adjusted R-squared: 0.07862
F-statistic: 2.294 on 6 and 85 DF, p-value: 0.04222
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: faith_pd
Df Sum Sq Mean Sq F value Pr(>F)
Village 4 127.26 31.814 1.0326 0.3964
Gender 1 59.86 59.857 1.9429 0.1676
Age_num 1 17.78 17.783 0.5772 0.4499
Residuals 72 2218.22 30.809
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-12.5271 -3.8100 -0.1492 3.3117 15.9465
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 25.03124 2.17450 11.511 <2e-16 ***
VillageKuyuwininha -2.67620 2.22460 -1.203 0.233
VillageShianana-Jiyakwanha -2.44842 2.32560 -1.053 0.296
VillageSudukuma -1.11259 2.42974 -0.458 0.648
VillageWashudihanha 1.66737 3.03047 0.550 0.584
GenderM 1.66382 1.29812 1.282 0.204
Age_num 0.02840 0.03738 0.760 0.450
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.551 on 72 degrees of freedom
Multiple R-squared: 0.08456, Adjusted R-squared: 0.008272
F-statistic: 1.108 on 6 and 72 DF, p-value: 0.366
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: faith_pd
Df Sum Sq Mean Sq F value Pr(>F)
Village 4 32.32 8.080 0.7172 0.5824
Gender 1 3.37 3.370 0.2991 0.5858
Age_num 1 222.48 222.483 19.7495 2.589e-05 ***
Residuals 87 980.08 11.265
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-6.2432 -1.7891 -0.2504 1.5857 10.1863
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.10518 1.54323 7.844 1.03e-11 ***
VillageKuyuwininha -1.41335 1.53533 -0.921 0.360
VillageMosenahanha 0.08616 1.69734 0.051 0.960
VillageShianana-Jiyakwanha -0.40834 1.53221 -0.267 0.790
VillageWashudihanha 0.32241 1.59235 0.202 0.840
GenderM -0.10671 0.71749 -0.149 0.882
Age_num 0.09539 0.02147 4.444 2.59e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.356 on 87 degrees of freedom
Multiple R-squared: 0.2085, Adjusted R-squared: 0.1539
F-statistic: 3.82 on 6 and 87 DF, p-value: 0.001995
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: faith_pd
Df Sum Sq Mean Sq F value Pr(>F)
Village 3 6.46 2.15 0.1136 0.95172
Gender 1 60.86 60.86 3.2110 0.07919 .
Age_num 1 507.04 507.04 26.7513 4.091e-06 ***
Residuals 50 947.69 18.95
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-9.407 -3.232 -0.059 2.255 11.458
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.12440 1.53300 6.604 2.50e-08 ***
VillageKuyuwininha -0.51548 1.61437 -0.319 0.751
VillageSudukuma 1.18120 1.83064 0.645 0.522
VillageWashudihanha -0.94686 1.97223 -0.480 0.633
GenderM 0.06533 1.36342 0.048 0.962
Age_num 0.19242 0.03720 5.172 4.09e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.354 on 50 degrees of freedom
Multiple R-squared: 0.3774, Adjusted R-squared: 0.3151
F-statistic: 6.061 on 5 and 50 DF, p-value: 0.0001848
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: faith_pd
Df Sum Sq Mean Sq F value Pr(>F)
Village 4 1439.7 359.92 3.7874 0.007236 **
Gender 1 20.3 20.33 0.2139 0.644995
Age_num 1 597.2 597.17 6.2839 0.014264 *
Residuals 78 7412.4 95.03
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-16.839 -6.341 -1.912 4.169 25.992
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.33590 4.63453 1.583 0.1175
VillageKuyuwininha 11.12761 4.78519 2.325 0.0226 *
VillageMosenahanha 1.96005 5.29088 0.370 0.7120
VillageShianana-Jiyakwanha 6.26461 4.82051 1.300 0.1976
VillageWashudihanha 10.13580 4.94157 2.051 0.0436 *
GenderM -2.13749 2.19252 -0.975 0.3326
Age_num 0.16052 0.06403 2.507 0.0143 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.748 on 78 degrees of freedom
Multiple R-squared: 0.2172, Adjusted R-squared: 0.157
F-statistic: 3.608 on 6 and 78 DF, p-value: 0.003288
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: faith_pd
Df Sum Sq Mean Sq F value Pr(>F)
Village 3 1000.7 333.57 3.3757 0.02678 *
Gender 1 14.0 13.95 0.1412 0.70895
Age_num 1 490.8 490.77 4.9665 0.03112 *
Residuals 43 4249.0 98.81
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-17.833 -6.404 -1.736 4.659 23.402
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15.9428 3.8103 4.184 0.000139 ***
VillageKuyuwininha 3.0606 4.0422 0.757 0.453084
VillageSudukuma -4.9909 4.4679 -1.117 0.270172
VillageWashudihanha 4.6992 5.7437 0.818 0.417782
GenderM -3.6879 3.3868 -1.089 0.282258
Age_num 0.2237 0.1004 2.229 0.031122 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.941 on 43 degrees of freedom
Multiple R-squared: 0.2616, Adjusted R-squared: 0.1758
F-statistic: 3.047 on 5 and 43 DF, p-value: 0.01928
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: faith_pd
Df Sum Sq Mean Sq F value Pr(>F)
Village 4 2227.4 556.86 5.6880 0.0004135 ***
Gender 1 388.1 388.08 3.9640 0.0496575 *
Age_num 1 40.9 40.94 0.4182 0.5195785
Residuals 86 8419.4 97.90
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-19.9594 -5.8154 -0.0413 4.7356 30.0663
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 54.70341 4.14874 13.186 < 2e-16 ***
VillageKuyuwininha -16.36338 4.24321 -3.856 0.000222 ***
VillageMosenahanha -17.91609 4.64464 -3.857 0.000221 ***
VillageShianana-Jiyakwanha -15.03259 4.22682 -3.556 0.000614 ***
VillageWashudihanha -10.90276 4.47693 -2.435 0.016944 *
GenderM 3.98777 2.11464 1.886 0.062699 .
Age_num 0.03924 0.06069 0.647 0.519579
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.894 on 86 degrees of freedom
Multiple R-squared: 0.2398, Adjusted R-squared: 0.1868
F-statistic: 4.522 on 6 and 86 DF, p-value: 0.0004973
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: faith_pd
Df Sum Sq Mean Sq F value Pr(>F)
Village 3 2057.4 685.81 3.8010 0.0220 *
Gender 1 22.3 22.28 0.1235 0.7281
Age_num 1 104.1 104.10 0.5770 0.4543
Residuals 26 4691.1 180.43
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-26.215 -8.914 3.249 8.816 31.260
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 52.7982 6.2340 8.469 5.97e-09 ***
VillageKuyuwininha -27.0667 9.3237 -2.903 0.00744 **
VillageSudukuma -5.3180 6.1586 -0.864 0.39575
VillageWashudihanha -17.0004 7.5285 -2.258 0.03256 *
GenderM -1.5252 5.9224 -0.258 0.79880
Age_num 0.1276 0.1680 0.760 0.45434
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 13.43 on 26 degrees of freedom
Multiple R-squared: 0.3176, Adjusted R-squared: 0.1864
F-statistic: 2.421 on 5 and 26 DF, p-value: 0.06271
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: faith_pd
Df Sum Sq Mean Sq F value Pr(>F)
Village 4 976.1 244.02 3.8654 0.006205 **
Gender 1 435.4 435.36 6.8964 0.010242 *
Age_num 1 7.7 7.70 0.1220 0.727727
Residuals 85 5366.0 63.13
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-14.9620 -4.7727 -0.1503 4.2880 24.7994
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 44.60976 3.49746 12.755 < 2e-16 ***
VillageKuyuwininha -9.63597 3.64022 -2.647 0.00967 **
VillageMosenahanha -8.83419 3.91219 -2.258 0.02650 *
VillageShianana-Jiyakwanha -4.03130 3.61637 -1.115 0.26811
VillageWashudihanha -6.79309 3.77159 -1.801 0.07523 .
GenderM 4.52324 1.70738 2.649 0.00962 **
Age_num -0.01731 0.04955 -0.349 0.72773
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.945 on 85 degrees of freedom
Multiple R-squared: 0.2092, Adjusted R-squared: 0.1533
F-statistic: 3.747 on 6 and 85 DF, p-value: 0.002348
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: faith_pd
Df Sum Sq Mean Sq F value Pr(>F)
Village 3 2271.4 757.12 5.7523 0.002796 **
Gender 1 101.5 101.49 0.7711 0.386227
Age_num 1 10.2 10.19 0.0774 0.782559
Residuals 33 4343.5 131.62
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-23.990 -5.218 -1.464 7.845 20.508
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 44.59878 3.99891 11.153 9.84e-13 ***
VillageKuyuwininha -22.14827 6.14464 -3.604 0.00102 **
VillageSudukuma 0.34254 4.40263 0.078 0.93845
VillageWashudihanha -5.75690 6.34527 -0.907 0.37084
GenderM 3.26711 4.06429 0.804 0.42723
Age_num 0.03458 0.12427 0.278 0.78256
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 11.47 on 33 degrees of freedom
Multiple R-squared: 0.3543, Adjusted R-squared: 0.2564
F-statistic: 3.621 on 5 and 33 DF, p-value: 0.01013
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: shannon
Df Sum Sq Mean Sq F value Pr(>F)
Village 4 2.993 0.74831 1.0878 0.3678
Gender 1 1.363 1.36264 1.9809 0.1629
Age_num 1 0.450 0.44994 0.6541 0.4209
Residuals 85 58.471 0.68789
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-2.0824 -0.6423 0.1359 0.6592 1.3825
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.870701 0.294116 19.961 <2e-16 ***
VillageKuyuwininha -0.174672 0.317054 -0.551 0.583
VillageMosenahanha 0.310353 0.360989 0.860 0.392
VillageShianana-Jiyakwanha -0.024449 0.308848 -0.079 0.937
VillageWashudihanha -0.215965 0.321496 -0.672 0.504
GenderM -0.277349 0.181881 -1.525 0.131
Age_num 0.004267 0.005276 0.809 0.421
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8294 on 85 degrees of freedom
Multiple R-squared: 0.07595, Adjusted R-squared: 0.01072
F-statistic: 1.164 on 6 and 85 DF, p-value: 0.3331
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: shannon
Df Sum Sq Mean Sq F value Pr(>F)
Village 4 3.922 0.98047 1.9663 0.1088
Gender 1 0.314 0.31397 0.6297 0.4301
Age_num 1 0.000 0.00001 0.0000 0.9968
Residuals 72 35.902 0.49864
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-1.84682 -0.42938 0.02862 0.48796 1.43152
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.460e+00 2.766e-01 19.736 <2e-16 ***
VillageKuyuwininha -2.237e-01 2.830e-01 -0.791 0.432
VillageShianana-Jiyakwanha 2.524e-01 2.959e-01 0.853 0.396
VillageSudukuma -1.078e-01 3.091e-01 -0.349 0.728
VillageWashudihanha 3.549e-01 3.855e-01 0.920 0.360
GenderM 1.298e-01 1.651e-01 0.786 0.434
Age_num 1.933e-05 4.755e-03 0.004 0.997
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7061 on 72 degrees of freedom
Multiple R-squared: 0.1055, Adjusted R-squared: 0.03099
F-statistic: 1.416 on 6 and 72 DF, p-value: 0.2204
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: shannon
Df Sum Sq Mean Sq F value Pr(>F)
Village 4 2.976 0.74390 1.1975 0.3177
Gender 1 0.512 0.51212 0.8244 0.3664
Age_num 1 1.091 1.09103 1.7564 0.1885
Residuals 87 54.044 0.62119
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-2.35162 -0.36635 0.02832 0.49628 1.73734
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.293606 0.362386 11.848 <2e-16 ***
VillageKuyuwininha -0.048434 0.360533 -0.134 0.893
VillageMosenahanha 0.477654 0.398575 1.198 0.234
VillageShianana-Jiyakwanha -0.031306 0.359799 -0.087 0.931
VillageWashudihanha 0.276176 0.373921 0.739 0.462
GenderM -0.185748 0.168483 -1.102 0.273
Age_num 0.006680 0.005041 1.325 0.189
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7882 on 87 degrees of freedom
Multiple R-squared: 0.07811, Adjusted R-squared: 0.01453
F-statistic: 1.228 on 6 and 87 DF, p-value: 0.2996
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: shannon
Df Sum Sq Mean Sq F value Pr(>F)
Village 3 2.998 0.9994 1.1125 0.3528618
Gender 1 1.019 1.0191 1.1344 0.2919507
Age_num 1 11.602 11.6019 12.9146 0.0007441 ***
Residuals 50 44.917 0.8983
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-1.97103 -0.65291 -0.08837 0.53085 2.28247
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.387193 0.333748 10.149 9.75e-14 ***
VillageKuyuwininha -0.271446 0.351462 -0.772 0.443551
VillageSudukuma 0.675170 0.398546 1.694 0.096471 .
VillageWashudihanha -0.254288 0.429371 -0.592 0.556362
GenderM -0.040747 0.296829 -0.137 0.891366
Age_num 0.029106 0.008099 3.594 0.000744 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9478 on 50 degrees of freedom
Multiple R-squared: 0.258, Adjusted R-squared: 0.1838
F-statistic: 3.477 on 5 and 50 DF, p-value: 0.008957
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: shannon
Df Sum Sq Mean Sq F value Pr(>F)
Village 4 15.040 3.7601 3.9173 0.005968 **
Gender 1 0.890 0.8903 0.9276 0.338471
Age_num 1 4.885 4.8855 5.0898 0.026868 *
Residuals 78 74.869 0.9599
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-2.6370 -0.5574 0.0151 0.5295 2.3070
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.556815 0.465776 5.489 4.85e-07 ***
VillageKuyuwininha 0.866755 0.480917 1.802 0.0754 .
VillageMosenahanha -0.213381 0.531740 -0.401 0.6893
VillageShianana-Jiyakwanha 0.318818 0.484468 0.658 0.5124
VillageWashudihanha 0.709996 0.496634 1.430 0.1568
GenderM -0.311185 0.220351 -1.412 0.1619
Age_num 0.014519 0.006436 2.256 0.0269 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9797 on 78 degrees of freedom
Multiple R-squared: 0.2175, Adjusted R-squared: 0.1574
F-statistic: 3.614 on 6 and 78 DF, p-value: 0.003246
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: shannon
Df Sum Sq Mean Sq F value Pr(>F)
Village 3 8.257 2.75239 2.5850 0.06547 .
Gender 1 0.587 0.58698 0.5513 0.46183
Age_num 1 2.009 2.00935 1.8872 0.17664
Residuals 43 45.784 1.06474
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-2.60663 -0.55941 0.05706 0.45811 2.48434
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.21546 0.39552 8.130 3.11e-10 ***
VillageKuyuwininha 0.25954 0.41959 0.619 0.539
VillageSudukuma -0.40686 0.46378 -0.877 0.385
VillageWashudihanha 0.81756 0.59622 1.371 0.177
GenderM -0.40555 0.35156 -1.154 0.255
Age_num 0.01431 0.01042 1.374 0.177
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.032 on 43 degrees of freedom
Multiple R-squared: 0.1916, Adjusted R-squared: 0.09763
F-statistic: 2.039 on 5 and 43 DF, p-value: 0.09217
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: shannon
Df Sum Sq Mean Sq F value Pr(>F)
Village 4 9.095 2.27378 5.0370 0.001075 **
Gender 1 2.925 2.92506 6.4797 0.012694 *
Age_num 1 0.015 0.01469 0.0325 0.857273
Residuals 86 38.822 0.45142
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-2.23705 -0.31392 0.06807 0.40145 1.70577
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.9816042 0.2817180 24.782 < 2e-16 ***
VillageKuyuwininha -0.6076504 0.2881329 -2.109 0.03786 *
VillageMosenahanha -1.0233288 0.3153916 -3.245 0.00168 **
VillageShianana-Jiyakwanha -0.8307218 0.2870199 -2.894 0.00481 **
VillageWashudihanha -0.2943695 0.3040035 -0.968 0.33561
GenderM 0.3587055 0.1435935 2.498 0.01439 *
Age_num 0.0007434 0.0041210 0.180 0.85727
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.6719 on 86 degrees of freedom
Multiple R-squared: 0.2366, Adjusted R-squared: 0.1834
F-statistic: 4.443 on 6 and 86 DF, p-value: 0.0005815
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: shannon
Df Sum Sq Mean Sq F value Pr(>F)
Village 3 4.9432 1.64773 1.8099 0.1701
Gender 1 0.8526 0.85265 0.9366 0.3421
Age_num 1 0.0988 0.09875 0.1085 0.7445
Residuals 26 23.6700 0.91038
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-1.9357 -0.3463 0.2187 0.5007 1.2221
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.947236 0.442820 15.689 8.97e-15 ***
VillageKuyuwininha -1.461267 0.662292 -2.206 0.0364 *
VillageSudukuma -0.043154 0.437464 -0.099 0.9222
VillageWashudihanha -0.357715 0.534773 -0.669 0.5094
GenderM -0.420837 0.420687 -1.000 0.3264
Age_num -0.003931 0.011935 -0.329 0.7445
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9541 on 26 degrees of freedom
Multiple R-squared: 0.1994, Adjusted R-squared: 0.04542
F-statistic: 1.295 on 5 and 26 DF, p-value: 0.2964
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: shannon
Df Sum Sq Mean Sq F value Pr(>F)
Village 4 3.490 0.87262 2.2542 0.06992 .
Gender 1 0.759 0.75883 1.9603 0.16512
Age_num 1 0.288 0.28827 0.7447 0.39059
Residuals 85 32.904 0.38711
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-1.76584 -0.36760 0.05787 0.40467 1.23096
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.356886 0.273876 23.211 <2e-16 ***
VillageKuyuwininha 0.033964 0.285056 0.119 0.9054
VillageMosenahanha -0.565078 0.306353 -1.845 0.0686 .
VillageShianana-Jiyakwanha -0.154536 0.283188 -0.546 0.5867
VillageWashudihanha -0.063571 0.295342 -0.215 0.8301
GenderM 0.200625 0.133700 1.501 0.1372
Age_num -0.003349 0.003880 -0.863 0.3906
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.6222 on 85 degrees of freedom
Multiple R-squared: 0.1212, Adjusted R-squared: 0.05916
F-statistic: 1.954 on 6 and 85 DF, p-value: 0.08139
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: shannon
Df Sum Sq Mean Sq F value Pr(>F)
Village 3 18.1479 6.0493 7.4850 0.0005919 ***
Gender 1 0.0752 0.0752 0.0931 0.7622392
Age_num 1 1.7271 1.7271 2.1370 0.1532384
Residuals 33 26.6704 0.8082
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-1.68378 -0.64642 0.01611 0.76874 1.56309
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.501204 0.313355 20.747 < 2e-16 ***
VillageKuyuwininha -1.641425 0.481495 -3.409 0.00174 **
VillageSudukuma 0.486595 0.344990 1.410 0.16776
VillageWashudihanha 0.033095 0.497216 0.067 0.94733
GenderM -0.001152 0.318479 -0.004 0.99713
Age_num -0.014235 0.009738 -1.462 0.15324
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.899 on 33 degrees of freedom
Multiple R-squared: 0.4279, Adjusted R-squared: 0.3412
F-statistic: 4.937 on 5 and 33 DF, p-value: 0.001753
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: observed_otus
Df Sum Sq Mean Sq F value Pr(>F)
Village 4 71933 17983.2 2.9929 0.02310 *
Gender 1 23138 23137.8 3.8508 0.05300 .
Age_num 1 16751 16750.6 2.7878 0.09867 .
Residuals 85 510733 6008.6
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-176.646 -55.980 6.271 53.262 178.467
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 281.7795 27.4882 10.251 <2e-16 ***
VillageKuyuwininha -51.6439 29.6320 -1.743 0.0850 .
VillageMosenahanha 7.4068 33.7382 0.220 0.8268
VillageShianana-Jiyakwanha -27.1231 28.8651 -0.940 0.3501
VillageWashudihanha -69.4344 30.0472 -2.311 0.0233 *
GenderM -37.7143 16.9987 -2.219 0.0292 *
Age_num 0.8233 0.4931 1.670 0.0987 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 77.52 on 85 degrees of freedom
Multiple R-squared: 0.1796, Adjusted R-squared: 0.1217
F-statistic: 3.102 on 6 and 85 DF, p-value: 0.008525
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: observed_otus
Df Sum Sq Mean Sq F value Pr(>F)
Village 4 19671 4917.7 1.4301 0.2328
Gender 1 9212 9212.3 2.6790 0.1060
Age_num 1 81 80.8 0.0235 0.8786
Residuals 72 247584 3438.7
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-121.01 -39.47 -2.29 38.95 196.28
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 213.77144 22.97306 9.305 5.66e-14 ***
VillageKuyuwininha -35.84516 23.50232 -1.525 0.132
VillageShianana-Jiyakwanha -14.89448 24.56937 -0.606 0.546
VillageSudukuma -17.21544 25.66953 -0.671 0.505
VillageWashudihanha 15.91064 32.01615 0.497 0.621
GenderM 21.97504 13.71427 1.602 0.113
Age_num 0.06054 0.39486 0.153 0.879
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 58.64 on 72 degrees of freedom
Multiple R-squared: 0.1047, Adjusted R-squared: 0.03013
F-statistic: 1.404 on 6 and 72 DF, p-value: 0.225
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: observed_otus
Df Sum Sq Mean Sq F value Pr(>F)
Village 4 5151 1287.7 1.0305 0.39621
Gender 1 135 134.6 0.1077 0.74358
Age_num 1 8307 8306.6 6.6474 0.01161 *
Residuals 87 108714 1249.6
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-57.02 -25.67 -2.26 20.37 137.47
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 90.6802 16.2534 5.579 2.69e-07 ***
VillageKuyuwininha -16.0785 16.1702 -0.994 0.3228
VillageMosenahanha 7.9291 17.8765 0.444 0.6585
VillageShianana-Jiyakwanha -10.3576 16.1373 -0.642 0.5227
VillageWashudihanha 3.5212 16.7707 0.210 0.8342
GenderM -5.4706 7.5566 -0.724 0.4710
Age_num 0.5829 0.2261 2.578 0.0116 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 35.35 on 87 degrees of freedom
Multiple R-squared: 0.1111, Adjusted R-squared: 0.04983
F-statistic: 1.813 on 6 and 87 DF, p-value: 0.1058
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: observed_otus
Df Sum Sq Mean Sq F value Pr(>F)
Village 3 905 301.5 0.1785 0.9104782
Gender 1 1151 1151.3 0.6816 0.4129697
Age_num 1 27029 27028.8 16.0015 0.0002091 ***
Residuals 50 84457 1689.1
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-81.170 -31.132 -1.022 17.458 100.691
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 68.4666 14.4720 4.731 1.87e-05 ***
VillageKuyuwininha -3.5870 15.2401 -0.235 0.814885
VillageSudukuma 19.0999 17.2818 1.105 0.274360
VillageWashudihanha -7.2194 18.6184 -0.388 0.699843
GenderM -6.3581 12.8711 -0.494 0.623483
Age_num 1.4049 0.3512 4.000 0.000209 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 41.1 on 50 degrees of freedom
Multiple R-squared: 0.2562, Adjusted R-squared: 0.1818
F-statistic: 3.444 on 5 and 50 DF, p-value: 0.009446
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: observed_otus
Df Sum Sq Mean Sq F value Pr(>F)
Village 4 75831 18957.7 3.8095 0.007003 **
Gender 1 2331 2330.6 0.4683 0.495789
Age_num 1 18054 18054.5 3.6280 0.060502 .
Residuals 78 388162 4976.4
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-111.15 -46.94 -12.22 21.41 216.82
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 31.5110 33.5376 0.940 0.3503
VillageKuyuwininha 81.4254 34.6279 2.351 0.0212 *
VillageMosenahanha 13.8680 38.2873 0.362 0.7182
VillageShianana-Jiyakwanha 39.0188 34.8835 1.119 0.2668
VillageWashudihanha 72.1992 35.7595 2.019 0.0469 *
GenderM -16.9188 15.8661 -1.066 0.2896
Age_num 0.8826 0.4634 1.905 0.0605 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 70.54 on 78 degrees of freedom
Multiple R-squared: 0.1986, Adjusted R-squared: 0.137
F-statistic: 3.222 on 6 and 78 DF, p-value: 0.007017
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: observed_otus
Df Sum Sq Mean Sq F value Pr(>F)
Village 3 49014 16337.9 3.6440 0.01988 *
Gender 1 5726 5725.9 1.2771 0.26471
Age_num 1 9324 9323.8 2.0796 0.15653
Residuals 43 192791 4483.5
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-113.62 -38.19 -14.26 30.16 213.71
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 88.7345 25.6661 3.457 0.00124 **
VillageKuyuwininha 27.0840 27.2280 0.995 0.32544
VillageSudukuma -34.4691 30.0953 -1.145 0.25841
VillageWashudihanha 54.5309 38.6894 1.409 0.16589
GenderM -35.1792 22.8132 -1.542 0.13039
Age_num 0.9749 0.6760 1.442 0.15653
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 66.96 on 43 degrees of freedom
Multiple R-squared: 0.2494, Adjusted R-squared: 0.1621
F-statistic: 2.858 on 5 and 43 DF, p-value: 0.02582
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: observed_otus
Df Sum Sq Mean Sq F value Pr(>F)
Village 4 138346 34587 4.6867 0.001806 **
Gender 1 33281 33281 4.5098 0.036573 *
Age_num 1 1294 1294 0.1753 0.676453
Residuals 86 634653 7380
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-193.647 -46.368 -6.585 40.051 299.562
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 386.9368 36.0200 10.742 < 2e-16 ***
VillageKuyuwininha -99.5728 36.8402 -2.703 0.008283 **
VillageMosenahanha -145.8888 40.3254 -3.618 0.000501 ***
VillageShianana-Jiyakwanha -86.3498 36.6979 -2.353 0.020905 *
VillageWashudihanha -63.6089 38.8694 -1.636 0.105394
GenderM 37.5990 18.3596 2.048 0.043618 *
Age_num 0.2206 0.5269 0.419 0.676453
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 85.91 on 86 degrees of freedom
Multiple R-squared: 0.2141, Adjusted R-squared: 0.1593
F-statistic: 3.905 on 6 and 86 DF, p-value: 0.001695
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: observed_otus
Df Sum Sq Mean Sq F value Pr(>F)
Village 3 203674 67891 4.3376 0.0132 *
Gender 1 4310 4310 0.2754 0.6042
Age_num 1 12620 12620 0.8063 0.3775
Residuals 26 406946 15652
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-222.88 -91.70 27.60 57.89 351.28
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 420.128 58.063 7.236 1.1e-07 ***
VillageKuyuwininha -282.647 86.840 -3.255 0.00314 **
VillageSudukuma -43.697 57.360 -0.762 0.45304
VillageWashudihanha -154.029 70.120 -2.197 0.03716 *
GenderM -22.778 55.160 -0.413 0.68303
Age_num 1.405 1.565 0.898 0.37746
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 125.1 on 26 degrees of freedom
Multiple R-squared: 0.3515, Adjusted R-squared: 0.2268
F-statistic: 2.819 on 5 and 26 DF, p-value: 0.03655
[1] "--------Below is ANOVA table for year 2015"
Analysis of Variance Table
Response: observed_otus
Df Sum Sq Mean Sq F value Pr(>F)
Village 4 51640 12910.1 2.9234 0.02565 *
Gender 1 15757 15757.1 3.5680 0.06231 .
Age_num 1 120 120.0 0.0272 0.86944
Residuals 85 375377 4416.2
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2015"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-127.519 -49.928 -0.568 39.532 214.100
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 305.41513 29.25251 10.441 <2e-16 ***
VillageKuyuwininha -38.85588 30.44654 -1.276 0.2054
VillageMosenahanha -66.76567 32.72127 -2.040 0.0444 *
VillageShianana-Jiyakwanha -4.86418 30.24703 -0.161 0.8726
VillageWashudihanha -29.78978 31.54525 -0.944 0.3477
GenderM 27.05143 14.28041 1.894 0.0616 .
Age_num -0.06833 0.41447 -0.165 0.8694
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 66.45 on 85 degrees of freedom
Multiple R-squared: 0.1524, Adjusted R-squared: 0.09262
F-statistic: 2.548 on 6 and 85 DF, p-value: 0.02564
[1] "--------Below is ANOVA table for year 2016"
Analysis of Variance Table
Response: observed_otus
Df Sum Sq Mean Sq F value Pr(>F)
Village 3 233404 77801 6.4756 0.001438 **
Gender 1 3103 3103 0.2582 0.614710
Age_num 1 365 365 0.0304 0.862612
Residuals 33 396478 12014
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] "--------Below is coefficient for year 2016"
Call:
lm(formula = paste0(AA, " ~ Village + Gender + Age_num"), data = Work_dat)
Residuals:
Min 1Q Median 3Q Max
-219.109 -56.237 -0.805 67.591 248.268
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 348.3303 38.2059 9.117 1.55e-10 ***
VillageKuyuwininha -207.8651 58.7065 -3.541 0.00121 **
VillageSudukuma 32.8143 42.0632 0.780 0.44088
VillageWashudihanha -41.4571 60.6233 -0.684 0.49885
GenderM 20.6938 38.8307 0.533 0.59766
Age_num -0.2071 1.1873 -0.174 0.86261
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 109.6 on 33 degrees of freedom
Multiple R-squared: 0.374, Adjusted R-squared: 0.2792
F-statistic: 3.943 on 5 and 33 DF, p-value: 0.006508lm_results_final <- do.call("rbind", lm_results)
anova_results_final <- do.call("rbind", anova_results)
# clip <- pipe('pbcopy', 'w') write.table(anova_results_final, file = clip,
# row.names = F, quote = F, sep = '\t') close(clip)